{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook uses `numpyro` and replicates experiments in references [1] which evaluates the performance of NUTS on various frameworks. The benchmark is run with CUDA 10.0 on a NVIDIA RTX 2070."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as onp\n",
    "from sklearn.datasets import fetch_covtype\n",
    "\n",
    "import jax.numpy as np\n",
    "from jax import random\n",
    "# NB: replace gpu by cpu to run this notebook in cpu\n",
    "from jax.config import config; config.update(\"jax_platform_name\", \"gpu\")\n",
    "from jax.tree_util import tree_map\n",
    "\n",
    "import numpyro.distributions as dist\n",
    "from numpyro.handlers import sample\n",
    "from numpyro.hmc_util import initialize_model\n",
    "from numpyro.mcmc import hmc\n",
    "from numpyro.util import fori_collect"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We do preprocessing steps as in [source code](https://github.com/google-research/google-research/blob/master/simple_probabilistic_programming/no_u_turn_sampler/logistic_regression.py) of reference [1]:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Data shape: (581012, 55)\n",
      "Label distribution: 211840 has label 1, 369172 has label 0\n"
     ]
    }
   ],
   "source": [
    "data = fetch_covtype()\n",
    "features = data.data\n",
    "labels = data.target\n",
    "\n",
    "# normalize features and add intercept\n",
    "features = (features - features.mean(0)) / features.std(0)\n",
    "features = np.hstack([features, np.ones((features.shape[0], 1))])\n",
    "\n",
    "# make binary feature\n",
    "_, counts = onp.unique(labels, return_counts=True)\n",
    "specific_category = np.argmax(counts)\n",
    "labels = (labels == specific_category)\n",
    "\n",
    "N, dim = features.shape\n",
    "print(\"Data shape:\", features.shape)\n",
    "print(\"Label distribution: {} has label 1, {} has label 0\"\n",
    "      .format(labels.sum(), N - labels.sum()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we construct the model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def model(data, labels):\n",
    "    coefs = sample('coefs', dist.norm(np.zeros(dim), np.ones(dim)))\n",
    "    logits = np.dot(data, coefs)\n",
    "    return sample('obs', dist.bernoulli(logits, is_logits=True), obs=labels)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Benchmark HMC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "time to compile sample_kernel: 2.6033902168273926\n"
     ]
    }
   ],
   "source": [
    "num_steps, num_samples = 10, 100\n",
    "step_size = np.sqrt(0.5 / N)\n",
    "trajectory_length = num_steps * step_size\n",
    "init_params = {\"coefs\": np.zeros(dim)}\n",
    "\n",
    "_, potential_fn, _ = initialize_model(random.PRNGKey(1), model, (features, labels,), {})\n",
    "init_kernel, sample_kernel = hmc(potential_fn, algo=\"HMC\")\n",
    "hmc_state, _, _ = init_kernel(init_params, num_warmup_steps=0, step_size=step_size,\n",
    "                              trajectory_length=trajectory_length,\n",
    "                              adapt_step_size=False, run_warmup=False)\n",
    "\n",
    "start = time.time()\n",
    "sample_kernel(hmc_state.update(step_size=1.))  # HACK: force fast compiling!\n",
    "print(\"time to compile sample_kernel:\", time.time() - start)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 100/100 [00:03<00:00, 27.10it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "number of leapfrog steps: 1000\n",
      "avg. time for each step : 0.0041951661109924316\n"
     ]
    }
   ],
   "source": [
    "start = time.time()\n",
    "hmc_states = fori_collect(num_samples, sample_kernel, hmc_state,\n",
    "                          transform=lambda state: {\"coefs\": state.z[\"coefs\"],\n",
    "                                                   \"num_steps\": state.num_steps},\n",
    "                          progbar=True)\n",
    "num_leapfrogs = np.sum(hmc_states[\"num_steps\"]).copy()\n",
    "print(\"number of leapfrog steps:\", num_leapfrogs)\n",
    "print(\"avg. time for each step :\", (time.time() - start) / num_leapfrogs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In CPU, we get `avg. time for each step : 0.041614791870117185`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD8CAYAAABw1c+bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJztnXl0JHd177+3qzf1on2bGY2k2Vd7PPZ4vO94A4ITQwAnwYRncPIChBBe8gLJCwdyXuLk8QgkEBwnGEwAQ56xMeQYY+OFMcaexbN59n2RZkZbS2qp9+X3/qj6VVfv1VJrq76fc3w86ip1V6m7v3Xre+/vXhJCgGEYhqkdbHN9AAzDMMzswsLPMAxTY7DwMwzD1Bgs/AzDMDUGCz/DMEyNwcLPMAxTY7DwMwzD1Bgs/AzDMDUGCz/DMEyNYZ/rAyhEa2ur6O3tnevDYBiGWTC89dZbw0KINjP7zkvh7+3txa5du+b6MBiGYRYMRHTW7L5s9TAMw9QYLPwMwzA1Bgs/wzBMjcHCzzAMU2Ow8DMMw9QYLPwMwzA1Bgs/wzBMjWEp4f+nl47jl8eG5vowGIZh5jWWEv5vvHoSvzrOws8wDFMKSwm/3UZIpnl4PMMwTCksJfyKQkix8DMMw5TEUsLPET/DMEx5LCX8io2QSrHwMwzDlMJSwm+32ZASLPwMwzClKCv8RLSUiF4hokNEdJCIPlVgHyKifyKiE0S0n4iuNGz7MBEd1/77cLVPwIhimxuP/8JYBFf/71/g+MDErL82wzBMpZiJ+JMAPiOEWA/gWgAfJ6L1OfvcC2CV9t/DAL4BAETUDODzAK4BsBXA54moqUrHnsdcefxHL01gaCKGI5dY+BmGmf+UFX4hxEUhxG7t3xMADgNYkrPbfQC+I1TeBNBIRIsA3A3gRSFEQAgxCuBFAPdU9QwMqBF/eqaevihDEzEAwHgkMeuvzTAMUykVefxE1AtgM4DtOZuWADhv+LlPe6zY44We+2Ei2kVEu4aGprYIS7ERknOQ3B2aZOFnGGbhYFr4icgH4EcA/kQIEaz2gQghHhNCbBFCbGlrMzU2Mg/7DNfxbz81gk1feAGjoXjW48Ms/AzDLCBMCT8ROaCK/veEEE8X2KUfwFLDz13aY8UenxEUm21GPf5TwyGMRxI4NRzKenx4Ur0QjIdZ+BmGmf+YqeohAN8EcFgI8eUiu/0EwINadc+1AMaFEBcB/BzAXUTUpCV179IemxHsM1zVE46nAACDwWjW48Ps8TMMs4Cwm9jnBgAfAvA2Ee3VHvscgG4AEEI8CuA5AO8EcAJAGMBHtG0BIvobADu13/uiECJQvcPPRrERkjOY3I0mVOEfyBV+zeoZi8TzfodhGGa+UVb4hRC/AkBl9hEAPl5k2+MAHp/S0VWI3UZIpGZO+CNaxH8pGMt6POPxJ6f9Gv/zqf1wO2z4wn0bp/1cDDPTjEcSCEYSWNrsmetDYSrAUit3lRmu448k8q2eRCqNUc3bD1bB6tl1NoC3+8en/TwMMxv844vH8M5/eg3BKNucCwlLCX+1PP7JWBInBvMXY0mPf2AiI/wBrcLH7bBhLDx9qycQimMyNv07B4aZDfpGI5iIJvEfb5yd60NhKsBSwq/YbFWp43/i12dw39deh8jp+5Px+DNWj1y8tbzVh1A8Zdpq2nNuNC9XkNTuHiaiLPzMwiAQUj//33r9tP79YOY/lhL+akX8Y+E4QvEUYslsEZcev1Gwpb+/st0HwJzdI4TA739rJ/755eNZj0vLaJKFn1kgBEJxLG2uw/BkHP+563z5X2DmBZYSfkWpTlWPFPxQjuUiPf6JaBLhuLpN1vBL4TdT0jkeSWA8ksCl8eyIf0SLnibjSaR5rgCzABgJxXHH2g5c1dOEf/3lqRktrmCqh6WEv1oRfyyhfnilpy+JGH6Wdo+M+Fe0qcI/ZkL4+0YjADI2kSSgXUSEAMJ828zMc+LJNCaiSTR7nfijW1egfyyCn+67MKXnSqbS+MovjuWtimdmBksJf7WqemJJVXQjOeIbSaTgsqt/Mmn3DE3EUOdQsKjRDcBcxF9M+IcNH/rp2j1j4Tge+dkRjsCYGUMWNjR7nbh9bTvWdvrx9VdOTMnr33t+DF/5xXG8eHig2ofJFMBSwl+1iD9ZJOJPpNDb4gWQEf7hyRha/U401DkAmPP4+0bDANTmbsYEcmAycyGYmGZ53EuHB/HoL0/i4IWqt1ViGAAZa7LF6wQR4X/esxYnh0L462cP5BVGlOP44CQAVKUyjimPpYS/Wr16dOHP9fjjKfS0qAtVBg1WT6vPpQv/mIl+PTLiT6RE1h1CwBDxT0yzpFOWnFZjbQHDFMIY8QPAbWvb8YnbVuI/d/XhyR2VJXqPD0xqz8mf19nAUsJfvYhfjfQLRfzt9S7UOZRMxD8RzxJ+M1ZP/1hE//ewIcqvptUjL0y8sIaZKaTwt/ic+mOfvnM1blndhs//5AB2nxs1/VzHtXUz7PHPDpYSfrUffxWqemRyN9fjj6fgcdrRUe/CwER2xO9QbPA6FdMev9epAAAGDT5/YDIOxaZ2x5juIq5BLeLnxnHMTDGiFSO0eF36Y4qN8NUPXoFFDXX4+Pd2m/4+ntCsngBbPbOCpYS/2h5/JJ4RXyEEIokU3A4F7fVuDIxHkUylEQjH0aZFPA11DpPCH8ampY0AshO8I6EYFmtJ4qpF/FXoH1Qp+86P4fvbz8366zKzSyCkBiryblfS6HHiD25ZjovjUX1IUSkmoglc1EqbOeKfHSwl/Godf/WsnlAsZXhMvRjUORR01LsxMBFFIByHEECbX4146uscZT3+8Yi6MndzdyHhj+vJ42p5/HMR8X9v+1l84acHK07wMQuLkVAcTR4HbLb8Ho4dfjWAya1cK4RM7HqcCkf8s4SlhL/qEb/B6pF+v8epoMPvwkAwiuEJ9UPa6lOFv6HOUTaZ2q8ldtctqodTsWVFRCOTcXRrXQ6nU9UjhJhTjz8QiiOWTLPNVIYXDw3gV8eH5/owpkwgFNMTu7nIYGgwWF74T2iJ3at6mkwVR1SbvtEw/vHFYzUVqFhK+GVVz3TfwMwCrkzULS8CdQ4FnQ1uRBNpnBpWP7Ct2oe80ZNv9Xxv+9msbp6ylHNpkwdtfpd+8UikVKFs87vgdSrTsnqCkaR+8aqG+D71Vh92nTE/RmFEu10fMPGlr2X+4fkj+D8vHJ3rw5gygVC8rPCbsXqOD07AZbfh8q4GjIXjMzpMqRDP7r2Ar750HOcDkfI7WwRLCb9du+Wc7uemUFWPXLXrdqoePwC9Rt4Y8RuHsfSPRfCXzxzA46+fyXoMALqa6tDqc+pfDOlttnid8Lnt00ruDhq6h1ajnPMfnj+Cb/36jOn9ZbXHpZwmdEw2A8Eojg9MLNj2HCOheFZi14is9DFr9axo86HV50JazH4J8rkRNRgzfm+sjqWEX1bETLdfT6aOPyP8UUPE36FFMxnhL5zcPavN5jVGy32jEbgdNjR7nWjzu/QvxoheGueCz2WflscvI+06h1KVL1Ewmqgo6RbQI/7a+SJVSjSRQjCaRDiewnntLnChUSrid9kVNHoc5oR/YBKrOnz6c822z382oH5PB00cq1WwlPDLiH+6t4q68Bfz+GXE3z8Ol90Gn0sdZNZQ50A0kdbvGM5okcT+vnH9wtE3GkZXkwdElC38k5nFMD63Y1pWjxTcle0+BKdZHRRLphBNpLMWl5VC9m8BgIHx6gp/Ki0wEIxiz7lRvHDw0oIebm8UxCOX8mc/zHeSqTTGwomiwg8A7X5X2Sg6FEuifyyCVe0+NHnU55rtyh5p8dRSoGJm5u6CIRPxT134k6m0fuGIFPD41XJONeIfCcWxpLEO6jx6oEH74I5HEmj3K3okEU+lsb9vHFuXNaN/LIKupjoAQJvPhUAohlRa6MvfW31O+F3TtXrU51rV7sMvjw1N+XkA6CJuVvhHDdFaNa2eSDyFW7/0Slbe4MHrevBFkyMqv/bycZwaDuHL77+iasc0HYzR5ZGLE7h7Q+ccHk3lyBbixsVbuRgDm2LI+v2V7f5MxD+Lwh9LpnBhXBV+jvgXKHrEP41hLMYe/IU8/jqHAo/TDr9bvWbKxC6AzOpd7UtxdjisJ7l2anZP32gESxo14fernuZIKGZY/u6C322fVlXPQDAKn8uOzgY3xiOJaSW7pVU0Go6beh5556Ieh7kvUjCawN/97HBeG2wju8+NYiAYwx/cshzf/PAW3LSqFT8/eMm0P/7mqQCe3t2PC2PzI4E3pEXCio1wdGDh9VPKbddQiDafq2xyV5ZyrurwodGjfn9GZ9Hq6R+NQH6szVQglUIIgc8+vR+vHh2swpHNLJYSfkVRT2c6EX/cIPyheAGPX1txK+2eNkPEk9u24WwgjMuWNGB1hw87zwQwGUtiLJxAV5NasimTwsMTcYxMxmEjoLHOAZ/LPi2rZ2gihvZ6F+rrHEimRV6X0UqQVlEiJUzdhcgvrcepmL513nZsCP/6y1P4r/3FW/ruOB0AEfDx21bijnUduP/KJRgIxrDf5HxiWdY61bbB1UZGl1csbcSRi7Nv9YyHE9NKKss71JLCr0X8pQKG44MTcCiEnmaPIeKfPQvvbEC1Y+02mnZy98ilCTy54zy+swDGUFpK+Kvh8RsjfqPVY/T4AaBDs3ukeAOqaAPQo+yzIyF0N3uwpbcZb50Z1asHdKvHUPI2oiXKbDaCzz3d5G4U7f7K+gcVw5gcHjXxhZRJ6nWL6k0Lv1zb8LMDl4rus/NMAOsX1aPerZ7T7Ws6YLcRfn6w+O8YkZbVs3vnh/APTcRgI+D6FS04MxKa1bGFQxMxXPfIS/jhNCZm6X16ilT1AEC7Xy17LhUwnBiYxPJWH+yKDXUOBS57dWZXm+W8JvwbljRMO+KXn8U3T43M+3bolhL+alT1yMSsjXKsHoPHD2QifqPwG4V2eDKOcDyF3hYPtvY2YyKWxEtar/E84Z+IYWQypn+JpMdvxlrZ3zeG//XjA1nR2+BEDB31bl0kp9O2wTj/10y1hWwtvW6RH8OTMVO9WmSJ6+snhgtepBKpNPacG8PVvc36Yw0eB65d3oIXTAp/MJKA22HDoYtBHB+Y+2TqYFDt8bR+UT3SItOdcjb4yb4LCMdTOGDybqkQpqweuYirhHd+fHASKzvUIUZEhGavc1Y9/nMjYbgdNmxYXD/tiP/nBwfgdtgQjqew9/xYlY5wZrCU8Fcz4m+oc2RN3DKWcwIGq6eAxz8WTuDsiJrY7Wn1YktvEwDgmb39AIAlmvDLi8bQRCyrNM7ntqtTuOLlo8AXDw3gP948i9Pa6wmhVr501LurE/Ebcg1ysHYpAqE4iIA1HX6kRWY0ZSn6tRLXRErg5SP5gzgO9I8jkkhh67LmrMfv2tCBk0MhPUFYDCEEJqJJvPvyxbBRftQfTaRmfdHQ4EQUbX4X1nT6AQCHL82ez//Mnj4AwNmRqZeRylxOk8dRdB9jYFOIiFbKurrdrz/W5HHOqsd/NhBGd7MHHX43RsMJPfCrlHMjYRy+GMTDNy2HjYDXprAi+9CFIN46G5iVdR2WEv5qVPXIVbtNHidCxqqeeAqKjeBQ1NeQtfzGiL/eILSylLOn2YOuJg8WN7hxaigEl92GNu13vC47PE5FjfhDcb1CwudSn2fChM8vRX1/nxphyFW77X4X6uvs2mPVsXrMeK+BcByNdQ4s1hLYZip7+sciuHFlKzrr3Xju7fwIXibGjRE/ANy5vgMA8MKh0lF/LJlGPJXG8jYvbljZimf39UMIgXgyjU//cC/W/q/nseJzz2HF557DHf/31VmxXQYnYmj3u9DT4oXbYcPRWSrpPD4wgQP9QTgVG85owcJUCITiaPQ4YFeKS0g54X+7fxxCZOZVA5j1iP98IIzuZq9u3ZpZd1AIafP89paluLyrEa+fqFz4H9t2Eg9/560pvX6lWEr47Tb1dKYX8atf+kaPWpMvr77heAoeh6KXbnY25Ef8io3gd9sxHkng3EgINoKeyL1ai1aXNGXKP+XvD09Kq0cVflkxNBkzN7gdAPadV2/b5e1qe4mI/1fHh023yzVG/Gbqq+Wdi7wjMuPz949G0NXkwT0bO7Ht2FBedc+O06NY1urN+lsDwKKGOmzqasALB0uP65Pn4Hc7cN8VS3A+EMHrJ0bw0BM78cyefnzo2h58+h2rcd8Vi3FyKFRRH/mpogq/G4qNsLrDXzXhHwxG8fB3dhUVsKf39EOxEd63pQsXxiJTjnBLLd6StPmKi6kQAl964SiaPA7cuKpVf7zJ69RLRWcaIQTOaRG/LNGeaknn8wcvYf2ieixt9uDGla3Ye36sosq8VFrgl8eGcMvqtoJN76qNpYRfj/irUM4pP9TS248kUnBriV0AuHVNO7543wZc1dOU9fty9e6ZkTAWN9bBqc3o3aJFq7KUU9Lmc+HCWATBaBIt2hfFpwn/VCJ+WULZ7ndlPH7DB/DopQn83je34xcmZ5sGI0k0ehxwKGTK4x+ZVJfxmxX+8UgCE7EkljTW4Z6NnYgl03jFUA6XTgvsOhvA1b1NBX//rg2d2Ht+DJdKLBaTOY56tx13b+iA027Df/v2Trx+Yhj/8L7L8Te/uRGfescqfOE9G6DYCL8+MVL2PAG1kd5D396J08OVRc6ptMDIZEwXmzUdfhypktXz+slhvHBoAP9ZIHGbTgs8u6cfN61qxZaeJqRFZhpcpQwbApViyM9NoZLO596+hB2nA/jMXWuy2jo3exyzFvHLPFx3cx3atW6iU0nwDk5EsfvcqL4W44aVrUilBbafMt/fau/5UYyGE7htbXvFrz8Vygo/ET1ORINEdKDI9iYieoaI9hPRDiLaaNj2aSI6SEQHiOhJInJX8+BzqY7HLyN+9UMtffZoIqX7+4Ca5H3wul79YiORwn82ENZbLAPAVk345R2ApM3vwlEt2SgvNn6XjPjLC7/sZnjwQhCJVFqP+Dvq3fqdgzHil1UMF02uqp2IJtBQ50CTx4mACb8+EIqjyetAi9cJxUZZwv/S4QHc/qVXs3InsqJnSVMdru5tRqvPmVXdc2JoEmPhRJ7NI7l7g2r3vFjC7pGRV73bAb/bgXs2dAIEPPp7V+H9W5bq+/ndDmzqasCvTN6m7zozipeODOK5ty+a2l8yEoohLdSLMwCsXVSP4cn4lG0GI9K3/9HuvrzigDdPj+DCeBS/tXkJerTP5pkKL1qSQIk+PRIiQpvPlSem0UQKf/vcYazt9OOBrd1Z25q8ToxHElUZqFSOc9oCy54WryHirzzB++KhAQgB3LNRFf4rexrhdthMf44A4JUjQ1BshJtXt1X8+lPBTMT/bQD3lNj+OQB7hRCXA3gQwFcBgIiWAPhjAFuEEBsBKAA+OK2jLYOiVKGqR/f41ShEduiMxLOFvxi68I+E9Pm8gLqK9qZVrbgl541t87v0yF72/JERv5la/mAkAafdhlgyjeMDk1kRv11R20kYq3ouGobE5/Ls3v68yDkYTcLvtqveq4mIfzQcR7PXBZuN0O534dJ45nV+fvASTg2HsqJbWdGzpLEOio1w5/pOvHJkUPfZt59Wo6bcxK5kRZsPazv9+NbrZ4qKhVyLIHMef3f/Zfjln92Kuwqslr1xZSv2942ZSojLqhh5t2UWKYRtWpS5VkvwVsPukcJ/aiiEfX3ZVTvP7O6Hz2XHXes70at9Ns9MMcEbCMXRXGLVrqTNn7+I69+2nUL/WAR//Rvr8wIn2bZhrMK81KO/PIkf7qxs+M85LQha2uxBi9cFG00t4n/+wCX0tniwWqtOctkVbF3WUpHP//KRQVzV05Q31GamKCv8QohtAErds6wH8LK27xEAvUTUoW2zA6gjIjsAD4AZLaKuZlVPbsQfTqT0xVulaPQ4cD4Qxlg4kSX8NhvhPx66Ro8KJMbkcLMWQcneP2Zq+ccjCf1uYn/fGAYn1FW7Xu056rWcg+SStjxdtoOWTEQT+NQP9uIHOV+eYCSBercDzV5nWY8/nRYYDSd0C6Cj3p0VQe05pwrkMUM5Zb/WoExWOr3rskUIx1P4xPf3YHAiip2nA2j3u/Q5BbkQET5z1xqcGg7hP3f1FdxHJqj9mvXlddmxqKGu4L43rGxFWqi12OWQTfr291VWFikje5mzkMJfDbvnzEgIG5fUw2W34UdvZf4eE9EEfnbgEu7d2Ik6p4JmrxN+t12vPqsE9X2Ol7V6gPy2DQPBKP7l1ZO4d2Mnrl/Rmrd/k7fyfj0nhybxD88fwQ93VrYu4dxIBERqebViI7T6yvcWymU0FMcbJ0dw94bOrNzdjStbcHxwsqAFGYols+7mL41HcehiELfPks0DVMfj3wfgfgAgoq0AegB0CSH6AXwJwDkAFwGMCyFeqMLrFaUqVT2a1SNtF93qqSDilwmiHoPVUwxjwlJW9fhNVvUIITAWSeCyrgbUu+3Y1zeOwWDGOwbUSiOjxy8tntyIX94p5D4ejKrC32Qi4h+PJJBKC/3L21Hv0j/445GEvjz/6KVM+WX/mFrKKUXkhpUt+Mt3rsO240O488vb8OrRQWxd1pz1pcrlHevasaWnCV/5xbGsGQoS+XeUOY9SbO5uQp1DMRWtHbgwDhupf9NKBENPwGvvfYvPhVafqyoR/7kRdbX4XRs68dP9F/TP8+efPYhwPIkPXdcDQL1g9rZ4pxTxj0USSIvSNfySNr87S/hfOHgJkUQKn7lrdcH9m2WjtgoSvF97+QTSovLE7NlACJ317qy1OcY2I6m0wC8ODZRcjPXT/ReQTAu854rFWY/fuFK9sy/0OfrD776F3/z667rlKXNaC034HwHQSER7AXwSwB4AKSJqAnAfgGUAFgPwEtHvFXsSInqYiHYR0a6hoak1FqtOVU+21SPfnIjJiL/ecKtmjPiL0WaI+Fu8lVk9obhaf97kceDyrkY94m/3Zwu/MeIfKGL1SDHKvRMIRpKor7Oj2VO+zE5vLa2dR2e9W389aYc4lOzeNP1jESw2NLojInzs5uV47o9vwoo2L4LRJK5d3lLydYkIf3HvWgxOxPAtw+wD/Rz0qp7yPQmddhuuWd5cVvjHwnH0jUZw2xr1y7r/vPmoP2P1ZN6ndYv8+h3EVJmIJjASiqO72Yv7r1yCsXACrxwZwo/39OPpPf345O2rcHlXo75/T4tnShF/wES7Bkmb34WRUGYh384zo+iod2FFm6/g/k1eh/Ya5iL+k0OTeHZvP9wOGwbLtIfI5XwgjKWGO0m1m2jme/HS4QF89Du78NfPHij6vD96qw9rO/3YsLgh6/G1nX60eJ15TRLHIwn8+uQITgxO4pGfHQag2jxLGuuwqr3w32QmmLbwCyGCQoiPCCGugOrxtwE4BeAdAE4LIYaEEAkATwO4vsTzPCaE2CKE2NLWNrUERzXr+KXVI2v5IwnzEb+kmD1hRH757TbSI1LFRvA4laxyzkMXgnjxUHYljhT0hjoHLutqwNFLEzgXCOsVNXKbsRY/E/Fnf7EycwGyLwgT0YzVUy7pJhfeSEFor3cjGE0iEk9hz7kxEAF3rO3IjvgNTeuMrGz34f/94fX47kPX4ANXL83bnsuW3mbcub4Dj756Mk80JqIJ/W9qhhtXtuLkUAgXx4tXvBzSRPr9Vy+FYqOKfP6hyRga6hx6pAkAV3Y34cil4LSa80l/v7fFg5tWtqLN78Jj207ir358AFf1NOGTt6/M2r+3xYu+0UjF7QXk4q1yyV1A/XwLkRHyt86OYktP8Ts4+dkxu4jrn186DpddwUduWIZ4Ml3RKvWzI2H0GIW/3qU3z5PHCgBP7jiPx7adyvv9E4MT2Nc3jvdd1ZW3zWYj3L62Ha8cHcz6+75+YhiptMDWZc144o2zePHQAF4/MYzb1raVvKutNtMWfiJqJCJ56f8ogG1CiCBUi+daIvKQekZ3ADg83dcrRcbjn37LBplk0iP+uEmPvy5jc3ic5SNM2d2zSevTI/HltGb+2ivH8dmn3876XdnTpKFOrUZJpgUGgrHsiN+dEX4hhG69DE9mR0e5cwEAtUV1KJ5CfZ0q/EKUXgVsnCkAqBE/oN5l7D43ilXtPmzpbdLXLQDIalOdi2Ij3LiqFY4Si4SM/PndaxCKJ/Fvr2V/SYMRNUFt9ot1w0rVe369RFnngQtqhH91bzNWtfuwtwKffzDnPQLU5HVaZMRmKkjh72nxwq7Y8JtXLMbuc2MgAF/5wBV5i616WjxIpYVeWWUWM+0aJPKOdnAihgtjEfSPRfSV7IWQ3zszEf+JwUn8ZN8FPHhdj54nMWu5ReIpDE7EsoKzNr8bI6G4LtS7z41i09JGvPvyRfi7nx3Bz3Kqt556S10Tcd8VSwq+xl0bOjERTWaVdb56dBB+tx2P//7VWNXuw8e/vxvheGpWbR7AXDnnkwDeALCGiPqI6CEi+kMi+kNtl3UADhDRUQD3AvgUAAghtgN4CsBuAG9rr/XYDJyDTrXq+IkyFSBho9VTQcTf01ze3wcylTy5iTKf257l8Z8ZDuu9+yWZiN+ZdQufF/FrzzMRUyc+tfpciCWzm2fJW1yjBSRf3++2Z5JuJSKxXEGQx3FxPIo958aweWmT3qLg6MAEookUhifjBSP+qbCqw4/LljTkWSbyrsUsazr8aPU5S9o9B/qDWNzgRrPXiU2azWbWZhiciGblYQBgc3cjFBvpq5SnglyJ261ZjB+4einqHAr+7r2XZVkakmWt6mf0dIV2T2ZanDmrB1DvcnZpF7UtPYUrtAC1TNrjVEwld7/+ygm47Aoevnl5pg7fpM8vp551G+zYjnr17mR4MoaENkPjqu4mfOm3N2FzdyP+5Id7sV1L+qfSAs/s6cMtq9vyFhZKblrVijqHoq8sF0JdpHXTqlb4XHZ85YNXQAgBl92G65bnJ7pnEjNVPQ8IIRYJIRxCiC4hxDeFEI8KIR7Vtr8hhFgthFgjhLhfCDFq+N3PCyHWCiE2CiE+JISY0UkHdqU6VT0uu02P1rPKOU1E/Lrwm/D3AbX0q6HOkfcl8rsywi87faZFtvAGDVbPoga3XiHUbhD++jr1ziGZSuvR/mVL6gFk2z2/M4yXAAAgAElEQVRyIHwwmtRbUwcN9e/NeiRWPOLP9X47G9Tj2X56BOORBDZ3N2JNhyr8xy5NZEo5i0T8U6HV59LvJiSyJNUsNhvh+hWt+NWJ4aJifvDCODYsUX3dy5c2YCycMD2se3AilpXbAQCP046Ni+ux87S5iP8XhwZw7d++lJW4PzcSRqs2uhNQh5sc+MLdePfliws+hyw+OFthLb+8wMvovBTthrYNu84E4HEqWLfIX/J3mjzlCwki8RSee/sifntLF1p8rorr8GWn3O4sjz+ziOvwxSBiybRWk6/g3x7cgiVNdfj9b6kL/14/MYyBYAzvvTLf5pG4HQpuXt2KFw4OQAiBwxcnMBCM4dbVanS/YXED/v69l+PP7l5jSluqiaVW7tqr4vGn4LIruh8ciacghNrT3m0i4pfDJMwKP6Am9lZ3ZH8ZjAPXhyfj+mwAoxUjF281eBwgImzqUoXIaCPIC9FENKn7+5dpgmWM7o2Rkvxi6yteNatH3Vb82h0IJeB1KvrfSV6AntcWZG3ubkKb34UmjwNHByYyi7cazf+tytHic2b9jYDKI35A9fmHJmI4VqBrZiiWxKnhEDYsVi+gm7S7rX0mfH4hhNquoT5/LePVvc3Y2zdmqo3Cm6dGcCkYxQ6DjXBmJKTX50ty6+SNtPqc8DqViit7AqE4/G67viq9FG1Zwj+Kzd2NJfv7ADBVOvzmqRHEkmm8Y51aOS4/8+Xq8IUQeOtsAN/foZYtGyvv9OeYiGG3dndyZbdqS7X6XPjhw9ehp8WDj3x7J/7++SOod9txx7rSFs1d6ztxKRjF2/3jeqL3ljWZHOb9V3bhozctL/kcM4GlhN9G1Yv4HYoNDoUQiqszZwGYSg4ubfKg1efC1mWlK1GMfPeha/BX71qf9ZjflZm7K1cYAtliLa0eOQdA2j25Hj+gRu9yBq6MVEdyhN+pfSHla2RWvNpNDckIhGJZi3r8WhO6I5cm4HPZsbLdB6JMb5qZiPhbfGoViTFSlx5/Jcj+Ma8dz68wO3wxCCGAjVolx5pOP5x2m6kEr7yjyvX4AbWfUzyZxtsm8gWyTcT205k8xNmRcJZ1UQ4iQk+Lt+LKniET7RokbocCv9uOk0OTOHIpWNLmkTR6HAiUKed89egg6hyKvrDP57KjzqGUtHp2nA7glv/zKt77jTfw+olhPLC1O6u7qLHNyJ7zY+isd+vNBgH1Ivbkx67F6g4fDl4I4jc2LS4bDN6+th2KjfDCwQG8enQQ6xbVZ1mxc4WlZu7Kcs7p1fGn4XKoz+Nx2hGJJ/V+PaY8fo8Du/7qHRW9ZqEIyBjxnxnORGS5wm83VKt8cOtS2Cjj3QLZMwJkxC8j1aEcq2eV9oGWHq5u9dQ5TI3FGwnFdUsIUIWlo96N08MhbFraoEefazr9eHp3P/pGw1BspHc6rQYtXicSKYFgJIkGj7zbSWSV2ZphcWMdVrb7sO34cF5EJnMIG7ULqEOxYf2i+ryVsoWQVSOFfOEtWt+nHWcCem+nYpzShP9NLeKPJlK4FIxmtQkxQ2+rp6IJYEII7D47qt/lmKHN78LLRwaRFvkdVgvR7HWWbBkthMArR4dw/YoWXXiJCO31rpLC/41XTyAcT+FLv70Jd2/o0Bf0SVp9ThBpEf859e4klyavE9/76LX4l1dO4MHre8ueS5PXia29zfjJvgu4MBbBx26e/ei+EJaK+BWlOlU9Lrv6YfI4FYTjqYqEv1r4XJm5u8aIzGhjjEfUPjqyWqWj3o1P3rEqq3pFCl4wksSlYAQtXic6690gAoa1L0k0kUIwmsTaznrtNWL67wBqctftUOB1KiWrLQp1bJTtbjcvzVRyrOn0YzKWxM7To+isd5e99a8EfZylwZKq1OOX3LSqFdtPjeS1aT7QP44Wr1M/NwDY1NWAA/3jSKUFdp0J4Pe/taPg3cKg3lIjP+pr8bmwos2LnadLJ3gTqTTOBdQBIgcvjKvdYAOyoqcy26y3xYvzo2HTvXFODE7i4ni0op4y7X4XxsIJ2Ai4ooCY5lKuJ//p4RDOBcK4dU32MbT7XXquKpd4Mo3tpwO4Z2MH3ndVV57oA2oA1uJ14tCFIM4HIrrNk0tDnQOffec600UJd67vwLlAGMm0wK2z1IunHJYS/ko9/ngynfeljmtWD6DO1w0nUnpJp3sWEzB+d2YK15mRMJY01sFuo6yIfyyS0KPaYhgj/kvjUXQ2qELb5HHqzyVLOWXSTV5cjBE/ADT7Si/iGg3F9bYTEnlbe2VP5gsvE7xvnRutqs0DZIRfnkMqrc4KrtTjB4CbV7Uhlkxj15nshOuBC0FsWNKQdYG9vKsR4XgKDz2xE+979A28enQIX3/lRN5zyog0t6pHsnVZM3adHS1pV54LhJFKC/zG5YuRFsCuM4GsUs5K6G3xIpESuDBmLikqfeqbV5uvQpE9idYvrtcTz6Vo9joxEU3qZZUjk7Gs4SSvHFWP4dY12f56e84qYSP7+sYQjqdwQ4E2EbnPIS/Yxs/sdJBzI/wuO67sKV7KOptYSviVCnv1/O1zh/Hhx3dkPRYzCL/HqSAcS+rC75nliD+tTeE6GwhjWasXzd7sxGVQi/hLoQ9jiapWz6IGOTIyI/xSjFa0+eC02/THg5EEiACfVuGUu3o3kUrrFUBCiKxhMhI5t+AKQ8S/ShP+VFqgq0qlnBL5+vKuZdJQklop1yxvhkOhrMg9lkzh+MAENmp2mWTTUlUkfn1iBB+/bQX+8JYVePNUQM9jSAZLWD2AaoVMRJNZ/YxyOT2k3gG+96ouOO02vHlqRL8rzE3ulqNHb9ZmzuffdnwYy9u8eV1mSyErmMz4+wCySocP9I/j+kdexv94ap++/dWjg1jR5s0rUW3zF7d6Xj8xDCLguhWlc2/t9Wqps0OhvNW4U2Vpszp+9Z6NnabXpMw08+MoqoS9wjr+/rGIfossiSXSBqvHnm31zGLEr7dtiCX1Tp+tWuJSMm5C+LMifm0kIyDLHlURN/rOrV6nXuYZjCbhd9n1hWXqkIyM8P/xk3vwoW9uB6BeoGLJdF6J34PX9eJrv7M5ywJqqHNgsXZBqHbEL4V/uECeolI8Tju29DRjm2GM3ltnR5FMizxRWNnuw9d+ZzNe+PTN+LO71+J3r1HbDf94T3/WfoPBGNwOm956OxfpgZeq5z81rFYareusxxVLG7H9dABnRkKod9v1Fedm6dXyQWYSvNFECttPjeDmVZXZFfIiV2rhlhGZJzofiOAT39+NtBB4enc/frynH+G4uiDqtjX51TTt9S5MxpIF+zX9+sQINi5uKPv3kUn39YsbTFXxmeX7H7sGj7z38qo933SxlPBXGvHHk+m8fjixZMqQ3FUQSaTyBq3PBtKD7BtVO332tnjR4nNmJWTHwgm9oqcYdQ4FdhthMBjDWDihR/wtPldexN9e79KrYgCtQVvWkIxMxB9NpPDykUFsPx3A233j+uO51R5LGusK1pGv1hZyVWvxlvEYAUOewrAWYSrctLoVhy8GMTgRRSKVxhd/egid9e6skjzJuy9frAvp0mYPru5twtM5ffGHJtXJW8VWEXc11aGz3o0dJXz+08MhtHid+sD5A/3jONAf1F+7Etr9LvjddhzoL98naMfpAGLJdF5r8XJIi+cak5Vusl/Pnz+1D+cCYTzx37ZiS08T/urHB/DDnecRT6XzbB71XAoPUwnHk9hzfhTXryz/+ro1aSIXUQl2xVaytHa2sZTwV1rVE0+mEYons76YeVZPPOPxz2ZyV0aEsoKkW0b8k5VF/ESEhjqHbh10NshB75nIfjAYg43U3ivGOni1DNIg/Ib66l1nRvWGdt9982xFy/iBjM9f7YhfzV849HPIdOacWgGbjG5/dXwYj207hSOXJvDF+zaY8qrvv7ILJ4dCeLs/U+1TqF2DESLCVb1N2Hu+eGnoyaGQXrl17XK11cPe82OmekMVer2bV7XhlaODZVcebzs2BKeiNrGrhFtWt2Hf5+8qam/lIj9DJ4dC+NM7V+P6Fa34ygevABHwxf86BI9TwdXL8u8ejHX4RnacDiCREmX9feNzFEvsWgVLCb9SYa+eRCqNtMiMVwSk8KsCX+ewqx5/QhUPs02+qoG0eg5qkVhvixetBlFOpwWC0fLCD6g2h5zy1WmweiZjSUQTKQxORNHqc0GxEVq8rqxo2SiYTV6ntq4hhddODMGhEN6zaTGe3deve8RmhnMA6mIuxUZFuzROh6y7lsjUrR4AWL+oHs1eJ36w4zy++tJx3Luxs+AAl0K887JFcNpteHq3avdcGo/ixNBk0cSuZEWrFxfGInr+JJfTwyEsb1OF/8ruJn39RaWlnJLb1rZjcCJWtjvotuNDuHpZk6keVLlUEu3K5m83rmzFf79VbSzX1eTB//6tyyAEcP2KVv07aqTY6t3XTwzDqdhMlZJu7m7C8lZv2VzAQsdSwl9xVY9WNWDsWaOu3FX/LF6XrOpR95tVj1+LKGUzsO5mD1p8LkQSKYRiSUxEkxACaDDh6dbXOfRqB31IvGEQtrqSVP251e/EcCgOIYQ6hKUuO+IH1KTba8eGcWV3E/7gluWIJtJ698Jmkx7z3Rs68Nqf35a1QKZatBjyFBPTSO4CavuGG1e2YseZAFx2G77wng2mf7ehzoE713Xgp/suYOeZAN7ztV8hHEvi967tKfl73S1ebR5ufi37RDSBoYkYlmsXTLdDwRVaYrnSUk7JrWvaQKS2By7GxfEIjg1MVuzvT4U2vwv/9uAWfP13rsy6YLxn02I8cv9l+NM7C/fyL2b1vH5iBJu7G019fzcuacDL/+PWrAFJVsRSwm+zEYgq8/gBIBTLifgdhnLO+Nx4/FL4jw1MoKPehTqnovvnI5PxrJbM5TBG7VL4W/3ac4XUWa/yQtDqdam5D+3iYvTGZeL2+MAkDl0M4qZVrdiwuAFX9TTp0aLZiJ+IZkT0gex+PdP1+AHo9eKfvXddwVYLpfitzUswEorj/f/6BlwOG57+oxsKTp4yIgX8bCBf+OWKXeMivWs166XSUk5Jq8+FTV2NeKmE8L92TE1wz9ZM2DvXdxQsVf7g1m6sz6mokjRpw91z248cuhjUO64yKpZauQuoUX8lHj+QPfDEaPV4HPasBPCsevyaWCdSQv9CyxbOw6EYHFo+w6zVA6h5A3lBkbfTw1rEL9sPZMoh41pyN/MRkdt+uk+doHmTFv196NoevHV2FA6FilarzCYtPqe++lhG/L4pRvyAGml21Ltx/RRu/29Z04buZg+6mz345wc266WKpZA94s8Oh4A12dtOaaWcK9oyIn//lV04MTSp92CaCnesbceXf3FMDQIKePG/PD6Edr9Lb388H9GHuxusnjdOqi0tbjCR2K0lLBXxA6qXaDriL2T1JFNZyV1A7UFjt5GpplTVwpg8lLXZrQaxHouowtZYZgEXkLk4yGgfyFxEBiaiGJnMWD0t0gKajOUtfJIR//MHL6GhzqG3LLj3sk40e51o9jpndZhEMVq86krRRCqNYCQBj1OZVv20XbHhhpWtUzo3h2LDS5+5Bd/96DWmRB9QrY46h1Iw4j81HIKNkFXD3tvqxb/87lXTsiJvW9sOIdQa+VySqTReOzaEm1fP7rCQqdBWn72I61cnhuF1KlltyxkLCr/dZpuC1aMKfyotkEiJTMTvUv8/HIrParQPqGIjX1NG/Ho0HqrU6skXfmkbHbs0gbQwzH/VHj8zHIIQ2d649PgnokncuLJV919ddgWfvXctHtjaPcWzrS7SxgqE4piYYruGalLpRUdtnubRWwcbOTU0iaXNnoLJzemwYXE9Oupd+vxXI3vPjyEYTea1SJiPqG0bVOEXQmDbsSHcsNL8MJ9aYe7vy6vMVCJ+OV5RXgiMdfwAEJiMz2q7BonPbUckkdI9X31x0kQMaa30zozwy30WGYRfdk08dFH15uWyepnUkk3AjMldtS8QIESme6Xkt7eUH484W+g21mRMHxa/0Ohu9ujvgZHTw6Esf79aEKmjAn+67yLiyXTW3e2rR4dgI+CmlQtD+Hdpi9+OD06ifyyCT+SMnGQsGfETkibLOXWPX4v4ZR90vVePQ70uBuYg4gcytfyyTM9lV8W64ohf8+k7cxKTbT4XDmudGaXVI6N62RbAKJqKjXS758Z5nCxrNeQp5kPEPxV6Wjw4Fwhn9agRQqilnK0zM5T7tjXtmIwldeGUvHpsEFd2N5XtCzUfaPe7MRpOIJ5M4xUtWb0Q7lRmG8sJfyURv2wCFdKFX4v4Dd05AXUA+WzW8EtkQtLYY71VW3E7Hk7AZbeZqjTKePzZVTQtPqd+0ZNWj9NuQ73brtflG5O7gFo5saw1v0/KfELmKUZCsbzVxwuF7hYv4sk0BgyJyoFgDOF4Csvaqh/xA+pdnNNuw3MHMrNlByeiONAfXDDiKQOYockYXj06hLWdfixqmJnqsYWM5YTfbiNTvXrSmp8PZKp6Ygkp/Jk6fkCN+GezlFPic6kDUIxRt2yuZmbVrkQOgDdaPepzuYr+W5YN5tokD9+8HJ+5q3Ad9XyhJS/iX3jCr1f2GHz+U0Nqj54VM2D1AGpvovs2LcZ/7urTx3Ru08o4C7VImI/IAObU0CR2ngksmOOebSwn/IpiLuKPG/qPT2p1/LrV48i2etJidks5JVt6m3HH2uwPrrqyNl6R8G9d1oy/ete6vFpmKfYNdY6sC1uLz6nf/eQK/weu7i46w3W+4HfZ4VRsGJ6Mq4vQFqjVAyArwXtS1vDPUMQPAH98xyoIIfC1V44DUKt8Wn0urF9UuHZ+viEXcT2zux/JtMBtC+ROZbZZeN+IMthtNlN1/EbhL2f1ALO7aldSaIVii8+JHWfiaA0nTJVyAqp9U2iupz6cPaduu8XQUz/X6lkIEBFatDujYDSxICP+xY11UGyEs4axm3vOjqLJ48jL1VSTpc0efODqpfjBjvP46I3L8drxYbxjXYfeoXW+I62e5w5chN89f/rfzzesF/Gb9PgThj4ok/Hs5K4zp44fmBvhL0Srz4XRcByBUNx0xF/0ubSyx9zeMcae+maakc1HWnxOXBiLIJESC/Li5VBs6Gqq062edFrgl8eGcMss1NJ/4rZVsNkIf/S93RiPJBaMvw+o5chEQDSRxs2r2riMswiW+6uYreopGPHnePweg+jNhdVTiFafE0KogzOmm7SUkX3uGECZHPW57FUdizibtHgzeYqFGPEDakmnFP63+8cxEorPimfd2eDGh67twaGLQbWMc9X8reDKRR2fqH5+F9IFa7ZZmN/qEpiN+I2dD/OtHunxGyL+eSP86oc6lkzrSdup0iYj/hyrR5ZDLkRvXNLic+rD5RfqefS0ePQBKa8eHQLR7PXK+e+3roDHqWBzd1PFw13mGvl5LjQzgVFZmN+IEpjt1WMUftnPJVPHr4q8YiO47DbEkul5Y/W0GKpvpm31aM+V25tFRkwLNVIGMt1Hgek1aJtLepq9CEaTGAvH8crRQWzqajQ972C6tPpc+NbvX73gRB8AVnX44HUpBQfaMyqWE37TEb9m9Xicir5yN5azcldujyXT8ybiN/rvDdP0rrubPfizu9fgNzZlV+noEf8C9MYlxr/TQj0PuX5jz/kx7Osbw5/cMbtltNcsX5iNzf7+vZfrK9uZwizMb0QJ7DabqTp+GfE3eZx6W+ZcqwdQa5tHw4l5E/Eb6+2nG40RET5+W/5ydnlXsVAjZSCnMmmBnocs6fzuG2chBHvWZpmLNTcLjZr3+I2rV3PLOYFMZc98ifjr3XY4FLWqY7pWTzEyEf/CFEwgO+JfqJaVHKX48tFBtHid02q7zDBGygo/ET1ORINEdKDI9iYieoaI9hPRDiLaaNjWSERPEdERIjpMRNdV8+ALYVcqq+pp9DgRT6aRSKURS2Qv4ALmn/ATkR7NzpQw17sdsNtowSZFgew7o4Vq9XicdrT7XRBCTVQulFp6Zv5jJuL/NoB7Smz/HIC9QojLATwI4KuGbV8F8LwQYi2ATQAOT/E4TWO6jl8T/mZtEVQolixo9UiLZ75YPUCm/t7sAq5KsdkIf/ObG/HBedJmeSrIiF+x0by5aE8Fafdw6wGmmpQVfiHENgCBErusB/Cytu8RAL1E1EFEDQBuBvBNbVtcCDE2/UMuTaVVPU2GHvNS+J2G2nWvNlh6PomHjPhnyuoBgAe2dmPdAlmmXwhZ/VLvts/74SGl6G72wkbAzQuolp6Z/1TjHngfgPsBvEZEWwH0AOgCkAIwBOBbRLQJwFsAPiWEyG8yXkXMRvxS5OVw8FA8qU/fMgrFfIz4ZTQ7k8K/0JEtrBeqvy/52M3LcMPKlgVZVsnMX6qR3H0EQCMR7QXwSQB7oIq+HcCVAL4hhNgMIATgL4o9CRE9TES7iGjX0NDQlA/GdK+enIg/FEsilkhn2TxAxuOfT5UCq9r96Gqq4+XoZWj1uRasvy9Z21mP+6/smuvDYCzGtL8VQogggI8AAKmh8mkApwB4APQJIbZruz6FEsIvhHgMwGMAsGXLlikX4Zr3+NV9pCUwGUupg9ZzBN6jWT1z0Y+/GB+7aRk+fH3PXB/GvKe72cMXR4YpwLSFn4gaAYSFEHEAHwWwTbsYBInoPBGtEUIcBXAHgEPTfb1ymO7Vo63SlROl1ORuqmjEP588frtiW7A9dGaTL79/04L29xlmpigr/ET0JIBbAbQSUR+AzwNwAIAQ4lEA6wA8QUQCwEEADxl+/ZMAvkdETqh3AR+p6tEXQLERUmYWcMmqHhnxa8ndosI/jyJ+xhzG9hYMw2QoK/xCiAfKbH8DQMG15EKIvQC2TO3QpoZax1/Jyl01+Tepe/zZAr+ooQ4ep7Ig57YyDMMUwnJqZr5Xj7pPY67V48iO+O+7YjFuWtWqe/0MwzALHcsZxZVU9TgVG5x29b/JeGGrx67Y0D6DE48YhmFmG8sJfyW9euSkLZ/Lrq/czbV6GIZhrIblhN/8BK6ULvxel4JQLIVYIr+qh2EYxmpYTuXMz9wVepdLr9OOiWgS8QJ1/AzDMFbDcsJvuldPKmP1+N1Gq8dyfxKGYZgsLKdyis0GIYB0GfGXyV0A8LrsWb16GIZhrIzlVM6u2Tflon414ldtHa/LXrSOn2EYxmpYTvgVbVhFOZ9fjfjVfX1Ozeox2D8MwzBWxXIqZ7fJiL90ZY+xnNPrsmNSJndZ+BmGsTiWUznTEX/KWMevIBTPH7vIMAxjRSyncpmIv7TwJ1KZ5K7P0IeHPX6GYayO5YRfsamnZMbjdxiqeiRs9TAMY3Usp3KFIv7nD1zCr08OZ+2X27JBwsLPMIzVsVzLSd3jN/Tk/8ovjmFRgxvXr8gMrI4Zk7uGzpu8cpdhGKtjufA2U8efqeqJJlJ68laSSGUqeNjqYRimlrCcyhWq6okl04jkCH88lfH4/W4WfoZhagfLqZxC+R5/NJFCOJ7M2i+3ZYOEq3oYhrE61hN+sxF/1gKujNhzHT/DMFbHcipXqFdPLJlGOJER/nRaIJkWXNXDMExNYjmVy9Txq8ndZCqNVFogbIj44yl1m/T46xwKtBsFtnoYhrE8lhN+vY5fK+eMJVWRjyfTuv0jhV9G90Sk+/wc8TMMY3Usp3K5Hn/UYPHIBG9cuxgYO3FKu4c9foZhrI7lVC535a6M+AHoCd6EFvHLqh4AhoifrR6GYayN5YQ/N+I3Cr/0+WXE7ygo/Jb7kzAMw2RhOZWza8ndTMSfsXpCJa0eNdJn4WcYxupYTuUyEb8q7rFEvtUTKyD8XqcdTrsNpC0AYxiGsSqWE/7cOv7s5G6Ox2+M+N12jvYZhqkJyiodET1ORINEdKDI9iYieoaI9hPRDiLamLNdIaI9RPRf1TroUlTi8RuTu2s6/FjR5puNQ2QYhplTzIS43wZwT4ntnwOwVwhxOYAHAXw1Z/unABye0tFNgWJ1/AAQSWgef4GI/w9uWYEff/yG2TpMhmGYOaOs8AshtgEIlNhlPYCXtX2PAOglog4AIKIuAO8C8O/TP1Rz5Ef8+VZPoYifYRimVqiG8u0DcD8AENFWAD0AurRtXwHw5wDShX+1+uRV9RRI7hby+BmGYWqFaijfIwAaiWgvgE8C2AMgRUTvBjAohHjLzJMQ0cNEtIuIdg0NDU35YHKreqIFIv5YgTp+hmGYWmHaoxeFEEEAHwEAUmshTwM4BeADAN5DRO8E4AZQT0TfFUL8XpHneQzAYwCwZcuW0pPSS5C3cjdRPLnLVTwMw9Qi01Y+ImokIqf240cBbBNCBIUQnxVCdAkhegF8EMDLxUS/mihK4aoep92GSLx4cpdhGKZWKBvxE9GTAG4F0EpEfQA+D8ABAEKIRwGsA/AEEQkABwE8NGNHa4L8Xj1qlN/kcWTq+Dm5yzBMDVNW+IUQD5TZ/gaA1WX2eRXAq5Uc2FTJ786pjlj0uuz6MBaO+BmGqWUsp3x2fRBLJuJ32W3wOBWEY9m9eji5yzBMLWI55ZOTtIxtmV0OGzwOe4HunNyXh2GY2sNywk9EsNsoq0mby67A41IQ0a0ewQ3ZGIapWSwn/IDq8xuTuy6HZvUYIn4X2zwMw9QollQ/u42QSmWSuy67gjqHXV+5G0+l4ODELsMwNYol1S8v4pfJXcMgFi7lZBimVrGk+tkVW9YCrozwy149gks5GYapWSypftkRfxouh4I6p4JYMo1UWqgRPws/wzA1iiXVL7uqJwW3FvEDQCSRQiyZ5hp+hmFqFkuqX+GIX12kHI4nEU9xxM8wTO1iSfVTI37ZnVNL7jq0iD+eQoLLORmGqWEsqX55Eb/B6gnFUhzxMwxT01hS/ew2m17HH0um4daSu4A6dzeeTHO7BoZhahZLCn/hOn7p8ae4qodhmJrGkupnV9SqnlRaIJESaq8eLeIPx1NIpNJw2pU5PkqGYQHslIsAAAriSURBVJi5wZLCLyN+OYTF5bBlrJ64Ws7JK3cZhqlVLKl+sqpHzts1JnfDcZncZY+fYZjaxJLCn4n4VeF3OxSDx5/kXj0Mw9Q0llQ/u03t1aNbPcaVu7rHb8lTZxiGKYsl1U9G/FHd6lHgUGxwKIRwgqt6GIapbSypfrJXjzHiB4A6h4JQLIlkWnCvHoZhahZLqp9iIyRT2R4/AHicdoxHEgDAET/DMDWLJdVPreM3VPU41NP0OBWMhTXh54ifYZgaxZLqp2jJ3Wgix+pxKhjTIn4XR/wMw9QollQ/e045p8surR4F4+E4ALDHzzBMzWJJ9VPkAq7c5K7Trkf87PEzDFOrWFL91Ig/nZ/cdSgIsvAzDFPjWFL99Ig/x+P3OBVoTTs5ucswTM1SVv2I6HEiGiSiA0W2NxHRM0S0n4h2ENFG7fGlRPQKER0iooNE9KlqH3wxpMcfTWZX9chGbQDg4IifYZgaxYz6fRvAPSW2fw7AXiHE5QAeBPBV7fEkgM8IIdYDuBbAx4lo/TSO1TSKNogllshP7kp49CLDMLVKWfUTQmwDECixy3oAL2v7HgHQS0QdQoiLQojd2uMTAA4DWDL9Qy6PXcm0ZXYoBMWmduKUA9cB9vgZhqldqqF++wDcDwBEtBVAD4Au4w5E1AtgM4DtVXi9smSqetJ6tA9kR/ws/AzD1CrVUL9HADQS0V4AnwSwB0BKbiQiH4AfAfgTIUSw2JMQ0cNEtIuIdg0NDU3rgDJVPamshVpeo8fPVg/DMDWKvfwupdHE/CMAQEQE4DSAU9rPDqii/z0hxNNlnucxAI8BwJYtW8R0jkmxEdICiMTTWcLPVg/DMEwVIn4iaiQip/bjRwFsE0IEtYvANwEcFkJ8ebqvUwkKqZ5+OJ6Ey1HE6uGIn2GYGqVsxE9ETwK4FUArEfUB+DwABwAIIR4FsA7AE0QkABwE8JD2qzcA+BCAtzUbCAA+J4R4rqpnUABFUYU/FE/lRPyGqh6O+BmGqVHKCr8Q4oEy298AsLrA478CMCeDbe1aFU84lhPxO9jjZxiGsaT6KTb1tHIjfg97/AzDMNYUfhnxh2LJolYPCz/DMLWKJdVPLtgKx5MF6/iJMhcHhmGYWsOSwp+J+FNwO4xWjyr8DsUGIhZ+hmFqE0sKv4z4I4lUVsQvrR7u08MwTC1jSQW0K5lo3mWI+J2KDYqN2N9nGKamsaQCyqoeILten4jgcSgs/AzD1DSWVEBj4tZtqN0HVLuHa/gZhqllpt2rZz6iGIQ/d4Wux6nAzsLPMEwNY0nht2cJf27Eb5+b5cQMwzDzBEsKf6mI3+tUkExPq/knwzDMgsaSwm83Jncd2cJ/+7p2JJIs/AzD1C6WFH5jxO/OsXr+6NaVs304DMMw8wpLZjmL1fEzDMMwFhV+pURyl2EYptaxpPDbSyR3GYZhah1LqqJSYgEXwzBMrWNJ4bcXadnAMAzDWFT4szx+Tu4yDMNkYUlVLLVyl2EYptaxpPCXWrnLMAxT61hSFY11/JzcZRiGycaSws8RP8MwTHEsqYpc1cMwDFMcS6qijPgVG3HvfYZhmBwsqYqyqsfN0T7DMEwellRGGfG7OLHLMAyThyWFX0b87O8zDMPkU1YZiehxIhokogNFtjcR0TNEtJ+IdhDRRsO2e4joKBGdIKK/qOaBl0Jh4WcYhimKGWX8NoB7Smz/HIC9QojLATwI4KsAQEQKgK8DuBfAegAPENH6aR2tSYgIio141S7DMEwBygq/EGIbgECJXdYDeFnb9wiAXiLqALAVwAkhxCkhRBzADwDcN/1DNodiI7i5Tw/DMEwe1VDGfQDuBwAi2gqgB0AXgCUAzhv269MeKwgRPUxEu4ho19DQ0LQPys4RP8MwTEGqIfyPAGgkor0APglgD4BUpU8ihHhMCLFFCLGlra1t2gel2Ig7czIMwxRg2sPWhRBBAB8BACIiAKcBnAJQB2CpYdcuAP3TfT2zqBE/Cz/DMEwu01ZGImokIqf240cBbNMuBjsBrCKiZdr2DwL4yXRfzyyKzcZ1/AzDMAUoG/ET0ZMAbgXQSkR9AD4PwAEAQohHAawD8AQRCQAHATykbUsS0ScA/ByAAuBxIcTBmTiJQthtBBe3a2AYhsmjrPALIR4os/0NAKuLbHsOwHNTO7Tp8ad3rsbyNu9cvDTDMMy8Ztoe/3zl/VcvLb8TwzBMDcJeCMMwTI3Bws8wDFNjsPAzDMPUGCz8DMMwNQYLP8MwTI3Bws8wDFNjsPAzDMPUGCz8DMMwNQYJIeb6GPIgoiEAZ6f4660Ahqt4OAuBWjxnoDbPuxbPGajN8670nHuEEKZaG89L4Z8ORLRLCLFlro9jNqnFcwZq87xr8ZyB2jzvmTxntnoYhmFqDBZ+hmGYGsOKwv/YXB/AHFCL5wzU5nnX4jkDtXneM3bOlvP4GYZhmNJYMeJnGIZhSmAZ4Seie4joKBGdIKK/mOvjmSmIaCkRvUJEh4joIBF9Snu8mYheJKLj2v+b5vpYqw0RKUS0h4j+S/t5GRFt197zHxpGgFoGbbTpU0R0hIgOE9F1Vn+viejT2mf7ABE9SURuK77XRPQ4EQ0S0QHDYwXfW1L5J+389xPRldN5bUsIPxEpAL4O4F4A6wE8QETr5/aoZowkgM8IIdYDuBbAx7Vz/QsALwkhVgF4SfvZanwKwGHDz38P4B+FECsBjEIb+2kxvgrgeSHEWgCboJ6/Zd9rIloC4I8BbBFCbIQ6tvWDsOZ7/W0A9+Q8Vuy9vRfAKu2/hwF8YzovbAnhB7AVwAkhxCkhRBzADwDcN8fHNCMIIS4KIXZr/56AKgRLoJ7vE9puTwD4zbk5wpmBiLoAvAvAv2s/E4DbATyl7WLFc24AcDOAbwKAECIuhBiDxd9rqJMB64jIDsAD4CIs+F4LIbYBCOQ8XOy9vQ/Ad4TKmwAaiWjRVF/bKsK/BMB5w8992mOWhoh6AWwGsB1AhxDiorbpEoCOOTqsmeIrAP4cQFr7uQXAmBAiqf1sxfd8GYAhAN/SLK5/JyIvLPxeCyH6AXwJwDmogj8O4C1Y/72WFHtvq6pxVhH+moOIfAB+BOBPhBBB4zahlmpZplyLiN4NYFAI8dZcH8ssYwdwJYBvCCE2Awghx9ax4HvdBDW6XQZgMQAv8u2QmmAm31urCH8/AON09S7tMUtCRA6oov89IcTT2sMD8tZP+//gXB3fDHADgPcQ0RmoNt7tUL3vRs0OAKz5nvcB6BNCbNd+fgrqhcDK7/U7AJwWQgwJIRIAnob6/lv9vZYUe2+rqnFWEf6dAFZpmX8n1GTQT+b4mGYEzdv+JoDDQogvGzb9BMCHtX9/GMCzs31sM4UQ4rNCiC4hRC/U9/ZlIcTvAngFwPu03Sx1zgAghLgE4DwRrdEeugPAIVj4vYZq8VxLRB7tsy7P2dLvtYFi7+1PADyoVfdcC2DcYAlVjhDCEv8BeCeAYwBOAvjLuT6eGTzPG6He/u0HsFf7751QPe+XABwH8AsAzXN9rDN0/rcC+C/t38sB7ABwAsD/A+Ca6+ObgfO9AsAu7f3+MYAmq7/XAL4A4AiAAwD+A4DLiu81gCeh5jESUO/uHir23gIgqJWLJwG8DbXqacqvzSt3GYZhagyrWD0MwzCMSVj4GYZhagwWfoZhmBqDhZ9hGKbGYOFnGIapMVj4GYZhagwWfoZhmBqDhZ9hGKbG+P9N2+mdVWbFOAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(hmc_states[\"coefs\"][:, 0]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Benchmark NUTS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_samples = 10\n",
    "step_size = 0.00167132\n",
    "init_params = {\"coefs\": np.array(\n",
    "    [+2.03420663e+00, -3.53567265e-02, -1.49223924e-01, -3.07049364e-01,\n",
    "     -1.00028366e-01, -1.46827862e-01, -1.64167881e-01, -4.20344204e-01,\n",
    "     +9.47479829e-02, -1.12681836e-02, +2.64442056e-01, -1.22087866e-01,\n",
    "     -6.00568838e-02, -3.79419506e-01, -1.06668741e-01, -2.97053963e-01,\n",
    "     -2.05253899e-01, -4.69537191e-02, -2.78072730e-02, -1.43250525e-01,\n",
    "     -6.77954629e-02, -4.34899796e-03, +5.90927452e-02, +7.23133609e-02,\n",
    "     +1.38526391e-02, -1.24497898e-01, -1.50733739e-02, -2.68872194e-02,\n",
    "     -1.80925727e-02, +3.47936489e-02, +4.03552800e-02, -9.98773426e-03,\n",
    "     +6.20188080e-02, +1.15002751e-01, +1.32145107e-01, +2.69109547e-01,\n",
    "     +2.45785132e-01, +1.19035013e-01, -2.59744357e-02, +9.94279515e-04,\n",
    "     +3.39266285e-02, -1.44057125e-02, -6.95222765e-02, -7.52013028e-02,\n",
    "     +1.21171586e-01, +2.29205526e-02, +1.47308692e-01, -8.34354162e-02,\n",
    "     -9.34122875e-02, -2.97472421e-02, -3.03937674e-01, -1.70958012e-01,\n",
    "     -1.59496680e-01, -1.88516974e-01, -1.20889175e+00])}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "time to compile sample_kernel: 7.5240113735198975\n"
     ]
    }
   ],
   "source": [
    "_, potential_fn, _ = initialize_model(random.PRNGKey(1), model, (features, labels,), {})\n",
    "init_kernel, sample_kernel = hmc(potential_fn, algo=\"NUTS\")\n",
    "hmc_state, _, _ = init_kernel(init_params, num_warmup_steps=0, step_size=step_size, \n",
    "                              adapt_step_size=False, run_warmup=False)\n",
    "\n",
    "start = time.time()\n",
    "sample_kernel(hmc_state.update(step_size=1.))  # HACK: force fast compiling!\n",
    "print(\"time to compile sample_kernel:\", time.time() - start)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 10/10 [00:12<00:00,  2.29s/it]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "number of leapfrog steps: 2726\n",
      "avg. time for each step : 0.004618817187370138\n"
     ]
    }
   ],
   "source": [
    "start = time.time()\n",
    "hmc_states = fori_collect(num_samples, sample_kernel, hmc_state,\n",
    "                          transform=lambda state: {\"coefs\": state.z[\"coefs\"],\n",
    "                                                   \"num_steps\": state.num_steps},\n",
    "                          progbar=True)\n",
    "num_leapfrogs = np.sum(hmc_states[\"num_steps\"]).copy()\n",
    "print(\"number of leapfrog steps:\", num_leapfrogs)\n",
    "print(\"avg. time for each step :\", (time.time() - start) / num_leapfrogs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In CPU, we get `avg. time for each step : 0.04056205542001939`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD8CAYAAABw1c+bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XlclXXax/HPxSaCInDADWTVXHJDccGtsk3T0pqasqZVx2xarKepaZrm6ZmZZ56pmakpW8fU9mmZFm0x07JSExcUVNxSQEUQZRFBENl+zx8cGyqRg5zDfTjner9evJRz/865L47y5ea6f/fvFmMMSimlvIeP1QUopZRqWxr8SinlZTT4lVLKy2jwK6WUl9HgV0opL6PBr5RSXkaDXymlvIwGv1JKeRkNfqWU8jJ+VhdwOhERESYuLs7qMpRSqt3YtGlTkTEm0pGxbhn8cXFxpKWlWV2GUkq1GyKy39Gx2upRSikvo8GvlFJeRoNfKaW8TLPBLyK9ROQrEdkhIttFZO5pxoiIzBORvSKyVUSG2R8fKiKp9udtFZFrXfFFKKWUcpwjJ3drgfuNMZtFpDOwSURWGGN2NBozGehj/xgFvGD/sxK4yRizR0R62p/7uTGm1LlfhlJKKUc1G/zGmEPAIfvfy0VkJxAFNA7+acBrpuGuLutEJFREehhjvmv0OvkicgSIBDT4lVLKIi3q8YtIHJAErP/Rpiggt9HnB+2PNX7uSCAAyGppkUoppZzH4eAXkU7A+8C9xpiyluxERHoArwO3GmPqmxgzW0TSRCStsLCwJS8PgDGGZ1fuITPvWIufq5RS3sSh4BcRfxpC/01jzAenGZIH9Gr0ebT9MUQkBPgU+J0xZl1T+zDGzDfGJBtjkiMjHbr47AeOnajhX+sPMPPVjRQcq2rx85VSyls4MqtHgIXATmPMk00M+wi4yT67ZzRwzBhzSEQCgA9p6P+/57SqTyM0KIAFN4/geFUts17bSGV1rSt3p5RS7ZYjR/xjgRuBiSKSYf+4TETmiMgc+5ilQDawF3gJ+JX98Z8DE4BbGj13qJO/hu8N6BnCvBlJ7Mgv4963M6ivN67alVJKtVvSMBHHvSQnJ5vWrNWzcE0Of/pkB3POS+Shyf2cWJlSSrknEdlkjEl2ZKxbLtLWWreNjSOr8DgvfpNFQmQwP0/u1fyTlFLKS3hk8IsIf7jiXHJLKnn4g230CgsiJdFmdVlKKeUWPHatHn9fH569fhhxEcHMeWMTOUUVVpeklFJuwWODH6BLR38W3TwCH4HbXtlIaWW11SUppZTlPDr4AWJsQcy/KZm8oyeY88YmqmtPe/2YUkp5DY8PfoARceE8fvUg1mWX8MjibbjjTCallGorHnly93SuTIomu7CCZ1buJTGyE7efl2h1SUopZQmvCX6A+y46h+zCCh5btou4iGAuPbe71SUppVSb84pWzyk+PsITPx/C4OhQ7n07Qxd0U0p5Ja8KfoBAf19eumk44cEBuqCbUsoreV3wA3TtHMiCm5M5XlXLzFd1QTellHfxyuAH6N8jhGeuT2LnIV3QTSnlXbw2+AEm9uvGI1MGsHzHYR7/fJfV5SilVJvwqlk9p3Pr2Diyi47zz2+ySYgI5toRMVaXpJRSLuX1wS8iPHr5uewvruR3H2bSKzyIMYkRVpellFIu49WtnlMaL+h2xxubyS48bnVJSinlMhr8dqcWdPP1EWa+mqYLuimlPJYGfyMxtiDm3zicvKMnuP11XdBNKeWZHLnZei8R+UpEdojIdhGZe5oxIiLzRGSviGwVkWGNti0TkVIR+cTZxbtCclw4f716MOtzSvjdh7qgm1LK8zhycrcWuN8Ys1lEOgObRGSFMWZHozGTgT72j1HAC/Y/Af4GBAG3O69s15qeFEV24XHmrdxLYtdOzNEF3ZRSHqTZI35jzCFjzGb738uBnUDUj4ZNA14zDdYBoSLSw/6cL4Fy55btevddfA5TB/fg8WW7WJZZYHU5SinlNC3q8YtIHJAErP/Rpiggt9HnB/npD4d2RUT4+zVDGBIdyr3vpLPtoC7oppTyDA4Hv4h0At4H7jXGlDm7EBGZLSJpIpJWWFjo7Jc/Kw0LuiVjC+7ArNd0QTellGdwKPhFxJ+G0H/TGPPBaYbkAb0afR5tf8xhxpj5xphkY0xyZGRkS57qUpGdO7DwlmQqTtYx89WNVJzUBd2UUu2bI7N6BFgI7DTGPNnEsI+Am+yze0YDx4wxh5xYp6X6dQ/hmRn2Bd3eyaBOF3RTSrVjjhzxjwVuBCaKSIb94zIRmSMic+xjlgLZwF7gJeBXp54sIquBfwMXishBEbnUuV9C27igX1d+P3UAK3Yc5q/LdEE3pVT71ex0TmPMGkCaGWOAO5vYNv7sSnM/t4yJI7uwgn+uyiY+IpjrRuqCbkqp9kev3G2BhgXdBjC+TwSPLM5k7d4iq0tSSqkW0+BvIT9fH567YRjxEcHMeWMTWbqgm1KqndHgPwshgf4sumUE/r4+zHxlI0crdEE3pVT7ocF/lnqFBzH/puHkH6tizhu6oJtSqv3Q4G+F4bHh/M2+oNvDuqCbUqqd8Po7cLXWtKFRZBVWMO/LPSRGduKO83VBN6WUe9Pgd4L7LupDTlEFjy/bRXxEEJMG9rC6JKWUapK2epxARPjb1YNJignl3ncy2Hqw1OqSlFKqSRr8ThLo78v8G+0Lur2aRk5RhdUlKaXUaWnwO1Fk5w4sumUEJ2vrmTpvNYvTW7ROnVJKtQkNfifr270zS+eOp3+PEO59J4P7392iK3oqpdyKBr8LRIV25O3Zo7lnYm8+SD/I5c+sITNPb+SilHIPGvwu4ufrw39d0pd/zRpNRXUtVz2/lkVrcnSuv1LKchr8LpaSaOOzuRMY3yeCP36yg1++lkaJLvGglLKQBn8bCA8OYMHNyTx6+QBWfVfE5KdXkZpVbHVZSikvpcHfRkSEW8fG88GvxhAc4Mf1C9bx5PLd1NbpGj9Kqbalwd/GBkZ14eO7x/GzYdHMW7mXGS+tI6/0hNVlKaW8iAa/BYI7+PH3a4bw1LVD2ZFfxmVPr2ZZZoHVZSmlvIQjN1vvJSJficgOEdkuInNPM0ZEZJ6I7BWRrSIyrNG2m0Vkj/3jZmd/Ae3Z9KQoPr1nPDHhQcx5YxO/X5xJVU2d1WUppTycI0f8tcD9xpgBwGjgThEZ8KMxk4E+9o/ZwAsAIhIOPAqMAkYCj4pImJNq9whxEcG8f8cYfjk+ntfX7Wf6c9+y90i51WUppTxYs8FvjDlkjNls/3s5sBOI+tGwacBrpsE6IFREegCXAiuMMSXGmKPACmCSU78CDxDg58Pvpgzg5VtHUFh+kqnPrOHtDQd0zr9SyiVa1OMXkTggCVj/o01RQG6jzw/aH2vq8dO99mwRSRORtMLCwpaU5TEu6NuVz+aOZ3hsGA99sI2730qnrKrG6rKUUh7G4eAXkU7A+8C9xpgyZxdijJlvjEk2xiRHRkY6++Xbja4hgbx+2ygeuLQvn2UWMGXeatIPHLW6LKWUB3Eo+EXEn4bQf9MY88FphuQBvRp9Hm1/rKnH1Rn4+Ah3XtCbd29Pob4ernkxlRe+zqK+Xls/SqnWc2RWjwALgZ3GmCebGPYRcJN9ds9o4Jgx5hDwOXCJiITZT+peYn9MOWB4bBhL547n0nO78/iyXdz88gaOlFdZXZZSqp1z5Ih/LHAjMFFEMuwfl4nIHBGZYx+zFMgG9gIvAb8CMMaUAH8CNto//mh/TDmoS0d/nr0+ib9cNYiN+0q47OnVfPOdd54DUUo5h7jjzJHk5GSTlpZmdRluZ8/hcu76Vzq7D5dz+4QE7r+kLwF+eg2eUgpEZJMxJtmRsZoa7Uifbp1ZctdYbhgVwz9XZXPNi2s5UFxpdVlKqXZGg7+dCfT35c9XDuKFG4aRU1TBZfNWsyRDz5crpRynwd9OTR7Ug6Vzx9O3e2fmvp3BA//eQmW13uJRKdU8Df52LDosiHdmj+buib15b/NBpj6zhu35eotHpdSZafC3c36+Ptx/SV/enDWKipO1XPncWl5du0+Xe1BKNUmD30OMSYxg6T3jGdcngkc/2s7s1zdxVG/xqJQ6DQ1+D2Lr1IGFNyfz+6kD+Hr3EabMW01ppYa/UuqHNPg9jIgwc1w8C28eQf6xKr3YSyn1Exr8Hmps7whCAv30pu5KqZ/Q4PdQvj7CyHgbqdka/EqpH9Lg92ApiTb2F1eSrzdzV0o1osHvwVISbADa7lFK/YAGvwfr170zYUH+rNXgV0o1osHvwXx8hFHxNtZlF+sFXUqp72nwe7gxvW3klZ4gt0T7/EqpBhr8Hu77Pn92kcWVKKXchQa/h+vdtRMRnTroCV6l1PccuefuIhE5IiKZTWwPE5EPRWSriGwQkYGNts0VkUwR2S4i9zqzcOUYEWF0Qjip2udXStk5csT/CjDpDNsfBjKMMYOBm4CnAew/AH4JjASGAFNFpHerqlVnJSXRxuGyk+QUVVhdilLKDTQb/MaYVcCZbpA+AFhpH7sLiBORbkB/YL0xptIYUwt8A1zV+pJVS/2nz6/tHqWUc3r8W7AHuoiMBGKBaCATGC8iNhEJAi4Dejlhf6qF4iOC6RaifX6lVANnBP9jQKiIZAB3A+lAnTFmJ/A4sBxYBmQAdU29iIjMFpE0EUkrLNQVJZ1JREhJsLEuu0T7/Eqp1ge/MabMGHOrMWYoDT3+SCDbvm2hMWa4MWYCcBT47gyvM98Yk2yMSY6MjGxtWepHUhJtFB0/yd4jx60uRSllsVYHv4iEikiA/dNZwCpjTJl9W1f7nzE0tIP+1dr9qbOTkhABoMs3KKXwa26AiLwFnA9EiMhB4FHAH8AY8yINJ3FfFREDbAdmNnr6+yJiA2qAO40xpc4tXzmqV3hHokI7kppVzM1j4qwuRylloWaD3xgzo5ntqcA5TWwbf5Z1KScTEVISbXyx8zD19QYfH7G6JKWURfTKXS+SkmCjtLKGXQXlVpeilLKQBr8XSUnU+fxKKQ1+r9IztCOxtiCdz6+Ul9Pg9zIpCTbW5xRTV6/z+ZXyVhr8XiYl0UZ5VS078susLkUpZRENfi+j6/MrpTT4vUzXkEASIoO1z6+UF9Pg90IpCTY25JRQU1dvdSlKKQto8HuhlEQbFdV1bMs7ZnUpSikLaPB7odGn+vza7lHKK2nwe6GITh3o260z6/RCLqW8kga/l0pJtJG27yjVtdrnV8rbaPB7qdEJNk7U1LHloC6YqpS30eD3UqMTwhHRPr9S3kiD30uFBgXQv3uIBr9SXkiD34ulJNrYdOAoVTVN3gpZKeWBNPi9WEqCjeraetIPaJ9fKW+iwe/FRiaE4yOQmuX56/YUHT/JvqIKq8tQyi00G/wiskhEjohIZhPbw0TkQxHZKiIbRGRgo233ich2EckUkbdEJNCZxavWCQn0Z2BUF6+4Mcs9b6VzwRNfc+/b6RworrS6HKUs5cgR/yvApDNsfxjIMMYMBm4CngYQkSjgHiDZGDMQ8AWua1W1yulSEmxk5JZyotpz+/z5pSdYm1XMkOhQPsss4MInv+bRJZkUlp+0ujSlLNFs8BtjVgElZxgyAFhpH7sLiBORbvZtfkBHEfEDgoD81pWrnG10oo2aOkPa/jP9E7dvH21p+G/39HVD+eaBC7gmuRdvrD/AeX/7iieX76a8qsbiCpWC2rr6Nvtt1Bk9/i3AVQAiMhKIBaKNMXnA34EDwCHgmDFmuRP2p5xoRFw4fj7i0dM6F6fnkRQTSqwtmO5dAvm/Kwex4r4JXNCvK/NW7mXCX79iwepsnd2kLJNfeoLrX1rPtfNTqThZ6/L9OSP4HwNCRSQDuBtIB+pEJAyYBsQDPYFgEflFUy8iIrNFJE1E0goLC51QlnJEpw5+DI723D7/7oJydhWUM21Izx88nhDZieeuH8ZHd41lYFQX/vfTnVz4xDf8Oy1Xb0up2tTy7QVMfno12/OP8eCkvgR38HP5Plsd/MaYMmPMrcaYoTT0+COBbOAiIMcYU2iMqQE+AMac4XXmG2OSjTHJkZGRrS1LtUBKoo2tB49xvA2ONNrakow8fH2EqT8K/lMGR4fy+sxRvDlrFLZOATzw3lYmPbWK5dsLMEZ/ACjXqaqp47+XZDL79U3EhAfxyT3juTIpuk323ergF5FQEQmwfzoLWGWMKaOhxTNaRIJERIALgZ2t3Z9yvpSECOrqDRv3eVafv77esCQjn3G9I4jo1OGMY8f2jmDJnWN5/oZh1NUbZr++iZ+9sJb1HvqbkLLW3iPlTH/uW15L3c+scfG8f8cY4iOC22z/jkznfAtIBfqKyEERmSkic0Rkjn1IfyBTRHYDk4G5AMaY9cB7wGZgm31f813wNahWGh4bhr+vsM7D+vybDxwlr/QE04ae/mj/x0SEywb1YPl9E/jLVYPIKz3BtfPXccvLG/Tm9MopjDG8uzGXy5/5liPlJ3n5lhE8MnUAAX5te0lVs80kY8yMZranAuc0se1R4NGzK021lY4BviT1CvO4Pv/ijDwC/X245NzuLXqen68PM0bGcGVSFK+s3cfzX+3lsnmrmTa0J/df3JcYW5CLKlaerKyqht99mMnHW/IZk2jjH9cOpVuINZc2uf4sgmoXRifaeHblHo6dqKFLR3+ry2m1mrp6Pt16iIsHdKfTWZ4sC/T3Zc55icwYEcM/V2Wx6NscPt16iOtHxXD3xD5Edj5z+0ipUzJyS7n7rc3kl1bxwKV9mXNeIr4+Ylk9umSDAhou5Ko3sCHHM/r8q/cUcrSyhukOtnnOpEuQPw9O6sc3D1zAtSN68ab9GoAnlu+mTK8BUGdQX2/45zdZXP3CWurr4d3bR3PnBb0tDX3Q4Fd2STGhBPj5eMx8/sXp+YQG+TO+j/NmiHULCeTPVw7ii/86j4n9uvLMyr2cp9cAqCYUlp/kllc28pfPdnHxgG4svWc8w2PDrS4L0OBXdoH+vgyP8Yw+f8XJWlbsOMyUQT1cctIsPiKYZ68fxsd3jfv+GoCJf/+ad9Nyqa3TW1k6aktuKb/7cBvLMg953LUTq/cUMvnp1azPLuZ/pw/k+RuG0SXIfVqo2uNX30tJtPHkiu84WlFNWHBA809wUyt2HOZETR3Tk6Jcup9B0V14feYo1u4t4vFlu3jwva3MX5XNA5f25ZIB3WiYxax+7GhFNX/9fDdvbzyAAG+uP0CcLYiZ4xO4elg0HQN8rS7xrNXU1fPE8u948Zss+nTtxJuzRtG3e2ery/oJPeJX3xuTaANgfU77PupfnJFHVGhHhseEtcn+xvSOYPGdY3nxF8OoN4bbX9/EVS+sZZ0H/PbkTPX1hrc3HGDiEw2/Hd02Np7031/Cc9cPo0tHf36/OJMxj33Jkyu+o+h4+1tAL7ekkmteTOXFb7KYMTKGj+4a55ahD3rErxoZHB1KR39fUrOKmTSwh9XlnJWi4ydZvaeI2RMS8GnDE2giwqSBPbiofzfe23SQp77Yw3Xz13HeOZE8OKkv5/bs0ma1uKPMvGM8sjiTjNxSRsaF88fp59KvewgAUwb34LJB3dm47yjzV2Ux78s9/PObLH42PJpZ4+JJiOxkcfXN+3hLPg9/sA0Enrt+GFMGu/f3jwa/+l6Anw/Jce27z790W0O/ePpQ17Z5muLn68N1I2OYnhTFa6n7eO6rLKbMW8MVQ3py/yXnEGtru6sz3cGxyhqeWLGbN9btJzw4gCd/PoQrk6J+0gYTEUbGhzMyPpy9R46zcE027206yFsbDnBR/27MnpBAcmyY27XPTlTX8YePt/P2xlySYkKZd10SvcLd/zoPccf1SJKTk01aWprVZXil57/ey1+X7SbtkYuaXebAHV31/LdUVtex7N4JVpcCwLETNcxflcXCNTnU1hlmjIzhrom9Lbtwp60YY3h/cx5/WbqTo5XV3Dg6lv+6pG+LrhEpLD/Ja6n7eH3dfkora0iKCWX2+AQuObe75dMhAXYVlHHXv9LJKjzOHeclct/F5+Dva133XEQ2GWOSHRqrwa8aSz9wlCufX8uz1ycxdXDr58C3pQPFlUz421f8ZlI/7jg/0epyfuBIWRXzVu7h7Q25iMDlQ3oyc1y8R7aAdh4q47+XZLJx31GSYkL507SBDIw6+6+zsrqW9zYdZMHqHA6UVBJrC2LWuHiuHt7LkhPBxhjeWH+AP32ygy4d/fnHz4cyrk9Em9fxYxr86qzV1tUz9I8rmDa0J3++cpDV5bTIsyv38Pfl3/HtQxOJCu1odTmndaC4kkXf5vBuWi6V1XWkJNiYOS6eif26tuk5CVcor6rhHyv28GrqPrp09OehSf24eni0076uunrD59sL+OeqbLbklhIW5M+No2O5MSWuza6iLq2s5jfvb+Xz7Yc575xInvj5ELf5zViDX7XKrS9vYH9xJSt/fb7VpTjMGMPF/1hFeFAA785JsbqcZh07UcM7Gw/wyrf7yD9WRUJEMLeOi+dnw6IICmhfp96MMXy0JZ///XQnRcdPcv3IGB64tC+hQa6ZEmyMIW3/UeavyuaLnYfx9/XhZ8OimTU+nkQXngjeuK+EuW+lc6T8JL+Z1I+Z4+Ld6od1S4K/ff0PU20iJdHGV7sLOVxW1W560dvzy9h75Dh/vnKg1aU4pEtHf2ZPSOTWsfF8llnAwtXZ/H5xJn//fDc3jIrhppQ4undx//d+z+Fyfr8kk3XZJQyO7sKCm5IZ0ivUpfsUEUbEhTMiLpyswuMsWJ3D+5t/eCJ4RJzzTgTX1Rue/2ov//jiO6LDgnj/jjEu/xpdTY/41U9sO3iMy59dw1PXDnX5RVDO8n9Ld/LytzlsePiidnnxmTGGTfuPsmB1Dst3FOAj8v15gNb0x12l4mQt877cw8I1OQR38OPBSX25bkSMZSddi46f5LXU/byeuo+jlTUM7RXK7AkJXNrKE8GHy6qY+3Y667JLuGJIT/585UA6B7rPFbiNaatHtUpdvSHpj8uZPLAHj1892OpymlVXbxj72EoGRoWw4OYRVpfTageKK3l5bQ7vbsylorqOUfHhzBqfwIVucB7AGMPSbQX86ZMdFJRVcW1yLx6c1Bebm/S5T1TX8d6mXBasyWF/cSUx4UHMHBfPNcnRLW6hrdx1mF//e2vDlM1p53LN8Gi3m07amAa/arVZr6bx3eFyVj14gdWlNGttVhHXv7SeZ2YkcXkTt1hsj8qqanhnQy6vrN1HXukJ4mxB3DYunquHtzzEnCGr8Dj/89F2Vu8pYkCPEP40fSDDY9vm6uiWqqs3rNjRcCI4/UApofYTwTc5cCL4ZG0dj3+2m0Xf5tCve2eevX4Yvbu6/0VkGvyq1RatyeGPn+xw6xkypzz0/lY+3pJP2iMXt+t1XppSW1fPsu0FLFidQ0ZuKSGBflw/Kpabx8TSo4vr/20qq2t5duVeXlqdTaC/L7++pC+/GB3rFnPpHZG2r4T5q7JZ8f2J4Chmjks4bZjnFFVw91ubycwr4+aUWH57WX8C/dvH/yk9uataLcW+bk9qVjFXD2+bG0CfjZO1dSzddohLz+3ukaEPDVcDTx3ck6mDe7Jp/1EWrslm/qosFqzOZsrgHswcF8/gaOefbDTGsHzHYf748Q7ySk9w1bAofju5f7u7AU1yXDjJceFkFx5n4Zoc+xXBuVzUvyu/HJ/AyPhwRIQPNh/k94sz8ffzYf6Nw1t857b2pNngF5FFwFTgiDHmJ1MmRCQMWAQkAlXAbcaYTBHpC7zTaGgC8N/GmKecUrlyqb7dOhMW5O/2wf/VrkLKqmqZ1k5OQrfW8NgwhscOJ7ekklfW7uOdjbksychnZFw4M8fHc1H/bk45Et9fXMGjH23n692F9O3WmXdvT2FkvHusJX+2EiI78ecrB3Hfxefweup+Xkvdxxc71zEkugtRYR1Zuq2AkXHhPHXdUHq6+W+5rdVsq0dEJgDHgdeaCP6/AceNMX8QkX7Ac8aYC380xhfIA0YZY/Y3V5S2etzDHW9sYuvBY6z5zQVue1Lrjjc2sXFfCet+eyF+Fl4ub5Xyqhre2ZjLy982nAeItQVx65g4rknuRfBZ3HKyqqaO57/O4sVvsgjw9eHei/pw85g4S5cicJUT1XW8v/kgC1Znc6Ckkrsn9uHuib3b7f8jp7Z6jDGrRCTuDEMGAI/Zx+4SkTgR6WaMOdxozIVAliOhr9xHSqKNzzILyC054ZY3GC+rquHLXUe4fmRMu/1mba3Ogf7MGp/ALWPi+Hz7YRauyeZ/Pt7Bkyu+Y8aoGG5OiXP46PXLnYf5n4+3k1tygiuG9OR3U/q3m+s4zkbHAF9+MTqWGSNjqKyuddtpmq7gjB7/FuAqYLWIjARigWigcfBfB7x1phcRkdnAbICYmBgnlKVaKyWhoc+/NquIGJv7/ZssyyyguraeaU64r2575+frw5TBPZgyuAeb9h9l0ZocXlqVzYLVOUwZ1HAeoKmLjnJLKvnDxzv4YudhenftxL9+OYoxidavPdNWfH3Eq0IfnBP8jwFPi0gGsA1IB76/AamIBABXAL8904sYY+YD86Gh1eOEulQr9e7aiYhOHUjNLua6ke4X/Esy8oi1BTG0nV9F6WwN5wHCyC2p5FX7eYCPtuQzIi6MmeMSuHhAw3mAqpo6XlqVzbNf7cXXR3hocj9uGxvvkttVKvfS6uA3xpQBtwJIQyM4B8huNGQysPlHrR/VDogIoxPCSc0qxhjjVn3+I2VVrM0q5u6JfdyqLnfSKzyIR6YOYO5FfXg37SAvf5vDnDc2ERMexFXDolicnse+4kouG9SdR6YM8PgTmuo/Wv2jXURC7Uf1ALOAVfYfBqfMoJk2j3JfKYk2jpSfJLuowupSfuCjLfkYg7Z5HNA50J+Z4+L5+tfn88INw4js3IGnvtiDiPDabSN5/obhGvpexpHpnG8B5wMRInIQeBTwBzDGvAj0B14VEQNsB2Y2em4wcDFwu9MrV23iVJ8/NavYpSsfttSSjHwGRXVxq5rcnZ+vD5MH9WDyoB7kllTSNaQDHfw889oHdWaOzOqZ0cz2VOCcJrZVALazK025g/iIYLqFNPT5fzE61upygIYSYaPwAAAMdUlEQVSlA7blHeORKf2tLqXdag+3B1Suo2dx1BmJCGMSI1if3dDndwdLMvK/v4uVUqrlNPhVs1ISbBQdr2bPkeNWl4IxhiUZeYxJtHn0HHOlXEmDXzWr8bo9Vtty8Bj7iyuZNtQ7lmhQyhU0+FWzeoUHERXa0S2Cf3F6HgF+Pkwa6LkLaCnlahr8yiEpiTbW5RRTX29dn7+2rp5PtuZzYb+uhHjZlZZKOZMGv3JISoKN0soadhaUNT/YRdZmFVN0vFrbPEq1kga/cog79PkXZ+TROdCP8/tGWlaDUp5Ag185pGdoR2JtQazLtib4T1TX8XlmAZcN7NFu7oiklLvS4FcOS0mwsT6nhDoL+vxf7jpMRXUd05J07r5SraXBrxyWkmijvKqW7fnH2nzfi9Pz6RbSgVHxeiG4Uq2lwa8c1njdnrZUWlnNN98d4YohPdvNDb6Vcmca/MphXUMCSYwMJrWN+/xLtxVQU2d0No9STqLBr1okJdHGxpwSaurq22yfizPySIwM5tyeIW22T6U8mQa/apGUhAgqquvYltc2ff680hNsyClh+tAoveGKUk6iwa9aZHRCONB2ff6Pt+QDaJtHKSfS4FctYuvUgb7dOrdZ8C9Oz2NYTCgxNl0/Xiln0eBXLZaSaCNtfwkna+tcup/dBeXsKijXo32lnEyDX7XY6AQbVTX1bMl1bZ9/cUYevj7ClME9XLofpbxNs8EvIotE5IiIZDaxPUxEPhSRrSKyQUQGNtoWKiLvicguEdkpIinOLF5ZY3RCOCKu7fPX1xs+yshnfJ8IIjp1cNl+lPJGjhzxvwJMOsP2h4EMY8xg4Cbg6UbbngaWGWP6AUOAnWdZp3IjoUEB9O8eQmp2kcv2senAUfJKTzBtqC7RoJSzNRv8xphVQMkZhgwAVtrH7gLiRKSbiHQBJgAL7duqjTGlrS9ZuYOURBubD5RSVeOaPv/i9Dw6+vtyyQC94YpSzuaMHv8W4CoAERkJxALRQDxQCLwsIukiskBEgpt6ERGZLSJpIpJWWFjohLKUK6Uk2KiurWfzgaNOf+3q2no+3XaIiwd0I7iDn9NfXylv54zgfwwIFZEM4G4gHagD/IBhwAvGmCSgAnioqRcxxsw3xiQbY5IjI3W9dXc3MiEcH4F1Lujzr95TSGlljbZ5lHKRVh9OGWPKgFsBpOHSyhwgGwgCDhpj1tuHvscZgl+1LyGB/gyK6uKSdXsWZ+QTFuTPhHP0AEApV2j1Eb995k6A/dNZwCpjTJkxpgDIFZG+9m0XAjtauz/lPkYn2sjILaWyutZpr1lxspYVOwqYMrgH/r4621gpV3BkOudbQCrQV0QOishMEZkjInPsQ/oDmSKyG5gMzG309LuBN0VkKzAU+D/nlq+slJJgo6bOkLbPeX3+5TsKqKqp14u2lHKhZls9xpgZzWxPBc5pYlsGkHx2pSl3NyIuHD8fITW72GltmcXp+USFdmR4TJhTXk8p9VP6u7Q6a8Ed/Bgc3cVpF3IVHT/Jmr1FTBvaEx+94YpSLqPBr1olJdHGtrxjHD/Z+j7/p1sPUVevN1xRytU0+FWrpCREUFdv2Jhzpmv8HLM4I49+3TvTt3tnJ1SmlGqKBr9qleGxYfj7SqundR4oriT9QCnTk/RoXylX0+BXrdIxwJekXmGt7vMvycgD4PIhetGWUq6mwa9abXSije35xzh2ouasnm+MYXFGHiPjw4kK7ejk6pRSP6bBr1ptTKKNegMbzrLPvz2/jKzCCqbrSV2l2oQGv2q1pJhQOvj5nHW7Z0lGHv6+wmWDdCVOpdqCBr9qtQ5+vgyPDWNtVsvX56+rN3y0JZ/zzulKaFBA809QSrWaBr9yipQEG7sKyimpqG7R89bnFHO47CTTk/SkrlJtRYNfOUVKog2A9S2c1rkkPZ/gAF8u6t/NFWUppU5Dg185xeDoUDr6+7ZoPn9VTR1LMw9x6cDuBPr7urA6pVRjGvzKKQL8fEiOa9l8/q93F1JeVauzeZRqYxr8ymlSEm3sOXKcwvKTDo1fkpFHRKcOjLG3iZRSbUODXzlNSkJDgK9zoN1TVlXDl7uOMHVwD/z0hitKtSn9jlNOMyiqC506+DnU51+2rYDq2npdm0cpC2jwK6fx8/VhRFyYQzdgX7IljzhbEEOiu7RBZUqpxhy59eIiETkiIplNbA8TkQ9FZKuIbBCRgY227RORbSKSISJpzixcuacxiRFkF1VwuKyqyTGHy6pYm1XMFUOjENEbrijV1hw54n8FmHSG7Q8DGcaYwcBNwNM/2n6BMWaoMUZvwegFTs3nP9Psno+35GMMTB+qF20pZYVmg98Yswo40+pbA4CV9rG7gDgR0atxvFT/HiGEBPqdcfmGJRn5DI7uQkJkpzasTCl1ijN6/FuAqwBEZCQQC0TbtxlguYhsEpHZTtiXcnO+PsKoBFuTJ3izCo+zLe8YV+i6+0pZxhnB/xgQKiIZwN1AOlBn3zbOGDMMmAzcKSITmnoREZktImkiklZYWOiEspRVUhJs5Jac4ODRyp9sW5Keh4+gwa+UhVod/MaYMmPMrcaYoTT0+COBbPu2PPufR4APgZFneJ35xphkY0xyZGRka8tSFmqqz2+MYcmWfMYkRtA1JNCK0pRSOCH4RSRURE6tpzsLWGWMKRORYBHpbB8TDFwCnHZmkPIsfbt1JizI/yftnozcUvYXV3KFntRVylJ+zQ0QkbeA84EIETkIPAr4AxhjXgT6A6+KiAG2AzPtT+0GfGifrucH/MsYs8zZX4ByPz4+wugEG+uyijHGfD9lc0lGPgF+PkwaqDdcUcpKzQa/MWZGM9tTgXNO83g2MOTsS1PtWUqijc8yCzhQUkmsLZjauno+2ZrPRf27EhLob3V5Snk1vXJXucSpdXtO9fm/zSqm6Hg1VwzRJRqUspoGv3KJ3l07EdGpw/d9/iXpeYQE+nFBPz1xr5TVNPiVS4gIKYk2UrOKOVFdx+fbC7hsUA86+OkNV5Symga/cpmUBBtHyk8yf1U2FdV1TNMbrijlFjT4lcucms//3Nd76R4SyKj4cIsrUkqBBr9yoThbEN1DAqmureeKoT3x8dGVOJVyBxr8ymVO9fkBpulFW0q5jWbn8SvVGjPHxRMfEcyAHiFWl6KUstPgVy41MKoLA6P0LltKuRNt9SillJfR4FdKKS+jwa+UUl5Gg18ppbyMBr9SSnkZDX6llPIyGvxKKeVlNPiVUsrLiDHG6hp+QkQKgf1n+fQIoMiJ5bRn+l78kL4fP6Tvx394wnsRa4xx6IYXbhn8rSEiacaYZKvrcAf6XvyQvh8/pO/Hf3jbe6GtHqWU8jIa/Eop5WU8MfjnW12AG9H34of0/fghfT/+w6veC4/r8SullDozTzziV0opdQYeE/wiMklEdovIXhF5yOp6rCQivUTkKxHZISLbRWSu1TVZTUR8RSRdRD6xuhariUioiLwnIrtEZKeIpFhdk5VE5D7790mmiLwlIoFW1+RqHhH8IuILPAdMBgYAM0RkgLVVWaoWuN8YMwAYDdzp5e8HwFxgp9VFuImngWXGmH7AELz4fRGRKOAeINkYMxDwBa6ztirX84jgB0YCe40x2caYauBtYJrFNVnGGHPIGLPZ/vdyGr6xo6ytyjoiEg1MARZYXYvVRKQLMAFYCGCMqTbGlFpbleX8gI4i4gcEAfkW1+NynhL8UUBuo88P4sVB15iIxAFJwHprK7HUU8CDQL3VhbiBeKAQeNne+logIsFWF2UVY0we8HfgAHAIOGaMWW5tVa7nKcGvTkNEOgHvA/caY8qsrscKIjIVOGKM2WR1LW7CDxgGvGCMSQIqAK89JyYiYTR0B+KBnkCwiPzC2qpcz1OCPw/o1ejzaPtjXktE/GkI/TeNMR9YXY+FxgJXiMg+GlqAE0XkDWtLstRB4KAx5tRvgO/R8IPAW10E5BhjCo0xNcAHwBiLa3I5Twn+jUAfEYkXkQAaTs58ZHFNlhERoaGHu9MY86TV9VjJGPNbY0y0MSaOhv8XK40xHn9E1xRjTAGQKyJ97Q9dCOywsCSrHQBGi0iQ/fvmQrzgZLef1QU4gzGmVkTuAj6n4az8ImPMdovLstJY4EZgm4hk2B972Biz1MKalPu4G3jTfpCUDdxqcT2WMcasF5H3gM00zIZLxwuu4tUrd5VSyst4SqtHKaWUgzT4lVLKy2jwK6WUl9HgV0opL6PBr5RSXkaDXymlvIwGv1JKeRkNfqWU8jL/D/SsEJ5FnBicAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(hmc_states[\"coefs\"][:, 0]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Comparing to HMC with similar trajectory length:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 10/10 [00:10<00:00,  1.06s/it]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "number of leapfrog steps: 3000\n",
      "avg. time for each step : 0.0035375065008799236\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "num_steps, num_samples = 300, 10\n",
    "step_size = np.sqrt(0.5 / N)\n",
    "trajectory_length = num_steps * step_size\n",
    "init_params = {\"coefs\": np.zeros(dim)}\n",
    "\n",
    "_, potential_fn, _ = initialize_model(random.PRNGKey(1), model, (features, labels,), {})\n",
    "init_kernel, sample_kernel = hmc(potential_fn, algo=\"HMC\")\n",
    "hmc_state, _, _ = init_kernel(init_params, num_warmup_steps=0, step_size=step_size,\n",
    "                              trajectory_length=trajectory_length,\n",
    "                              adapt_step_size=False, run_warmup=False)\n",
    "\n",
    "sample_kernel(hmc_state.update(step_size=1.))\n",
    "\n",
    "start = time.time()\n",
    "hmc_states = fori_collect(num_samples, sample_kernel, hmc_state,\n",
    "                          transform=lambda state: {\"coefs\": state.z[\"coefs\"],\n",
    "                                                   \"num_steps\": state.num_steps},\n",
    "                          progbar=True)\n",
    "num_leapfrogs = np.sum(hmc_states[\"num_steps\"]).copy()\n",
    "print(\"number of leapfrog steps:\", num_leapfrogs)\n",
    "print(\"avg. time for each step :\", (time.time() - start) / num_leapfrogs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In CPU, we get `avg. time for each step : 0.04001419194539388`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Average time for each leapfrog (verlet) step"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "|               | HMC (n=100, l=10) | HMC (n=10, l=300) | NUTS (n=10, l~300) |\n",
    "| ------------- |:-----------------:|:-----------------:|:------------------:|\n",
    "| Edward2 (CPU) |                   |                   |      68.4 ms       |\n",
    "| Edward2 (GPU) |                   |                   |       9.7 ms       |\n",
    "| Numpyro (CPU) |      41.6 ms      |      40.0 ms      |      40.6 ms       |\n",
    "| Numpyro (GPU) |      4.2 ms      |       3.5 ms      |       4.6 ms       |\n",
    "\n",
    "*Note:* Edward 2 number is obtained from reference [1]."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Some takeaways:**\n",
    "+ The overhead of iterative NUTS is small. So most of computation time is indeed spent for evaluating potential function and its gradient.\n",
    "+ The average speed of long trajectory is better than short trajectory.\n",
    "+ GPU outperforms CPU by a large margin. The data is large, so evaluating potential function in GPU is clearly faster than doing so in CPU.\n",
    "+ Iterative NUTS is 1.7x faster in CPU and 2x faster in GPU than the reported speed in reference [2]. This illustates the win of a graph-mode (built using iterative algorithm) NUTS over an eager-mode (built using recursive algorithm) NUTS."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### References\n",
    "\n",
    "1. `Simple, Distributed, and Accelerated Probabilistic Programming,` [arxiv](https://arxiv.org/abs/1811.02091)<br/>\n",
    "Dustin Tran, Matthew D. Hoffman, Dave Moore, Christopher Suter, Srinivas Vasudevan, Alexey Radul, Matthew Johnson, Rif A. Saurous"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Environment (conda_numpyro)",
   "language": "python",
   "name": "conda_numpyro"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
